{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple Oscillator Example\n",
    "\n",
    "This example shows the most simple way of using a solver.\n",
    "We solve free vibration of a simple oscillator:\n",
    "$$m \\ddot{u} + k u = 0,\\quad u(0) = u_0,\\quad \\dot{u}(0) = \\dot{u}_0$$\n",
    "using the CVODE solver. An analytical solution exists, given by\n",
    "$$u(t) = u_0 \\cos\\left(\\sqrt{\\frac{k}{m}} t\\right)+\\frac{\\dot{u}_0}{\\sqrt{\\frac{k}{m}}} \\sin\\left(\\sqrt{\\frac{k}{m}} t\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scikits.odes import ode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#data of the oscillator\n",
    "k = 4.0\n",
    "m = 1.0\n",
    "#initial position and speed data on t=0, x[0] = u, x[1] = \\dot{u}, xp = \\dot{x}\n",
    "initx = [1, 0.1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need a first order system, so convert the second order system \n",
    "$$m \\ddot{u} + k u = 0,\\quad u(0) = u_0,\\quad \\dot{u}(0) = \\dot{u}_0$$\n",
    "into \n",
    "$$\\left\\{ \\begin{array}{l}\n",
    "  \\dot u = v\\\\\n",
    "  \\dot v = \\ddot u = -\\frac{ku}{m}\n",
    "\\end{array} \\right.$$\n",
    "You need to define a function that computes the right hand side of above equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rhseqn(t, x, xdot):\n",
    "    \"\"\" we create rhs equations for the problem\"\"\"\n",
    "    xdot[0] = x[1]\n",
    "    xdot[1] = - k/m * x[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve the ODE you define an ode object, specify the solver to use, here cvode, and pass the right hand side function. You request the solution at specific timepoints by passing an array of times to the solve member."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "solver = ode('cvode', rhseqn, old_api=False)\n",
    "solution = solver.solve([0., 1., 2.], initx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "   t        Solution          Exact\n",
      "------------------------------------\n",
      "   0               1               1\n",
      "   1       -0.370694       -0.370682\n",
      "   2       -0.691508       -0.691484\n"
     ]
    }
   ],
   "source": [
    "print('\\n   t        Solution          Exact')\n",
    "print('------------------------------------')\n",
    "for t, u in zip(solution.values.t, solution.values.y):\n",
    "    print('{0:>4.0f} {1:15.6g} {2:15.6g}'.format(t, u[0], \n",
    "           initx[0]*np.cos(np.sqrt(k/m)*t)+initx[1]*np.sin(np.sqrt(k/m)*t)/np.sqrt(k/m)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can continue the solver by passing further times. Calling the solve routine reinits the solver, so you can restart at whatever time. To continue from the last computed solution, pass the last obtained time and solution. \n",
    "\n",
    "**Note:** The solver performes better if it can take into account history information, so avoid calling solve to continue computation!\n",
    "\n",
    "In general, you must check for errors using the errors output of solve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error:  Could not reach endpoint Error at time 24.5780834078\n",
      "Computed Solutions:\n",
      "\n",
      "   t        Solution          Exact\n",
      "------------------------------------\n",
      "   2       -0.691508       -0.691484\n"
     ]
    }
   ],
   "source": [
    "#Solve over the next hour by continuation\n",
    "times = np.linspace(0, 3600, 61)\n",
    "times[0] = solution.values.t[-1]\n",
    "solution = solver.solve(times, solution.values.y[-1])\n",
    "if solution.errors.t:\n",
    "    print ('Error: ', solution.message, 'Error at time', solution.errors.t)\n",
    "print ('Computed Solutions:')\n",
    "print('\\n   t        Solution          Exact')\n",
    "print('------------------------------------')\n",
    "for t, u in zip(solution.values.t, solution.values.y):\n",
    "    print('{0:>4.0f} {1:15.6g} {2:15.6g}'.format(t, u[0], \n",
    "           initx[0]*np.cos(np.sqrt(k/m)*t)+initx[1]*np.sin(np.sqrt(k/m)*t)/np.sqrt(k/m)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution fails at a time around 24 seconds. Erros can be due to many things. Here however the reason is simple: we try to make too large jumps in time output. Increasing the allowed steps the solver can take will fix this. This is the **max_steps** option of cvode:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed Solutions:\n",
      "\n",
      "   t        Solution          Exact\n",
      "------------------------------------\n",
      "   2       -0.691508       -0.691484\n",
      "  60        0.843074        0.843212\n",
      " 120        0.372884        0.373054\n",
      " 180       -0.235749       -0.235745\n",
      " 240       -0.756553       -0.756932\n",
      " 300       -0.996027       -0.996814\n",
      " 360       -0.865262       -0.866242\n",
      " 420       -0.412897       -0.413742\n",
      " 480        0.192583        0.192521\n",
      " 540        0.726263        0.727236\n",
      " 600        0.989879        0.991682\n",
      " 660        0.885441        0.887581\n",
      " 720        0.452113        0.453622\n",
      " 780       -0.149122       -0.148921\n",
      " 840       -0.694753       -0.696119\n",
      " 900       -0.981996       -0.984613\n",
      " 960       -0.904317       -0.907187\n",
      "1020       -0.490547       -0.492616\n",
      "1080        0.105433         0.10503\n",
      "1140        0.661996        0.663643\n",
      "1200        0.972376         0.97562\n",
      "1260        0.921229        0.925021\n",
      "1320        0.527762        0.530648\n",
      "1380      -0.0616772      -0.0609338\n",
      "1440       -0.627855        -0.62987\n",
      "1500       -0.960354       -0.964723\n",
      "1560       -0.935939       -0.941048\n",
      "1620       -0.563703       -0.567643\n",
      "1680       0.0179025       0.0167187\n",
      "1740        0.592593        0.594867\n",
      "1800        0.946886        0.951941\n",
      "1860        0.949087        0.955237\n",
      "1920         0.59865         0.60353\n",
      "1980        0.025868       0.0275291\n",
      "2040       -0.556352       -0.558703\n",
      "2100       -0.931643         -0.9373\n",
      "2160       -0.960605        -0.96756\n",
      "2220       -0.632579       -0.638239\n",
      "2280      -0.0695761      -0.0717232\n",
      "2340         0.51911        0.521447\n",
      "2400        0.914706        0.920828\n",
      "2460         0.97026        0.977994\n",
      "2520         0.66523          0.6717\n",
      "2580        0.113083        0.115777\n",
      "2640        -0.48092       -0.483173\n",
      "2700       -0.896022       -0.902558\n",
      "2760       -0.978027       -0.986518\n",
      "2820        -0.69649        -0.70385\n",
      "2880       -0.156248       -0.159605\n",
      "2940        0.441822        0.443956\n",
      "3000        0.875509        0.882526\n",
      "3060        0.983777        0.993115\n",
      "3120         0.72638        0.734626\n",
      "3180        0.199071        0.203121\n",
      "3240       -0.401994       -0.403871\n",
      "3300       -0.853534       -0.860769\n",
      "3360       -0.987807       -0.997773\n",
      "3420        -0.75483       -0.763966\n",
      "3480       -0.241495       -0.246241\n",
      "3540        0.361435        0.362997\n",
      "3600         0.82981        0.837332\n"
     ]
    }
   ],
   "source": [
    "solver = ode('cvode', rhseqn, old_api=False, max_steps=5000)\n",
    "solution = solver.solve(times, solution.values.y[-1])\n",
    "if solution.errors.t:\n",
    "    print ('Error: ', solution.message, 'Error at time', solution.errors.t)\n",
    "print ('Computed Solutions:')\n",
    "print('\\n   t        Solution          Exact')\n",
    "print('------------------------------------')\n",
    "for t, u in zip(solution.values.t, solution.values.y):\n",
    "    print('{0:>4.0f} {1:15.6g} {2:15.6g}'.format(t, u[0], \n",
    "           initx[0]*np.cos(np.sqrt(k/m)*t)+initx[1]*np.sin(np.sqrt(k/m)*t)/np.sqrt(k/m)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To plot the simple oscillator, we show a (t,x) plot of the solution. Doing this over 60 seconds can be done as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEPCAYAAABRHfM8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXmcZVdxJvhFZlVWZW2qfa9SaQVJSLbAFmJ1whgQpg3D\nIozswRhvuHvoZjzuaWyGtkr2THu67Zn2YGyZcdsyw9iYxYYRBgmEUYIAt1gkpFJJhSS01ZqVtWRV\nZmUtWZnRf5x3lDdv3uUscZf3Kr7fL3+Z772b8eIsEV9EnHPPJWaGQqFQKBRl6GtaAYVCoVB0B5Qw\nFAqFQuEEJQyFQqFQOEEJQ6FQKBROUMJQKBQKhROUMBQKhULhhEYJg4j+iohGiGhXzudDRHSCiB7s\n/Hy4bh0VCoVCYbCg4e+/A8CfAPh/C675OjO/uSZ9FAqFQpGDRjMMZr4PwPGSy6gOXRQKhUJRjLav\nYTCAlxPRQ0T0JSK6ummFFAqF4kJF0yWpMjwAYBszTxLRGwF8HsCVDeukUCgUFyRaTRjMPJ74+y4i\n+jMiWs3Mx5LXEZEeiKVQKBQBYGbnsn+rS1JEtIGIqPP3DQAoTRYWzNyzP7feemvjOmjbtH3avt77\n8UWjGQYRfRLATwFYS0R7AdwKYCEAMPPHALwDwL8kovMAJgG8qyldFQqF4kJHo4TBzLeUfP6nAP60\nJnUUCoVCUYBWl6QUBkNDQ02rUBl6uW2Atq/b0evt8wWF1LHaBiLiXmiHQqFQ1AkiAvfKordCoVAo\n2gMlDIVCoVA4QQlDoVAoFE5QwlAoFAqFE5QwFAqFQuEEJQyFQqFQOEEJQ6FQKBROUMJQKBQKhROU\nMBQKhULhBCUMhUKhUDhBCUOhUCgUTlDCUCgUCoUTlDAUCoVC4QQlDIVCoVA4QQlDoVAoFE5QwlAo\nFAqFE5QwFAqFQuEEJQyFQqFQOEEJQ6FQKBROUMJQKBQKhROUMBQKhULhBCUMhUKhUDhBCUOhUCgU\nTlDCUCgUCoUTlDAUCoVC4QQlDIVCoVA4oVHCIKK/IqIRItpVcM1HiOgJInqIiK6vUz+FQqFQzKLp\nDOMOADflfUhEPwPgcma+AsCvA7i9LsUUCoVCMReNEgYz3wfgeMElbwbw8c619wNYSUQb6tBNoVAo\nFHPRdIZRhi0A9iZe7wOwtSFdFAqF4oJG2wkDACj1mqW/YGwM+MM/BM6ejZPzF38B3H13nIxDh4A/\n+iPg/Pk4Od/4BvCpT8XJmJwEfuVXgHPn4uTcfTfwrW/FyWAGfvM3gZGRODmf+xxwxx1xMmZmgI9/\nHNi3L07O/fcDH/5wnAxm4IMfBI4di5Ozdy9wu0DB93d/F3jssTgZDz0EfOIT8bp85jPAN78ZJ+PZ\nZ4F/9+/idfnbvwV++MM4GadOAe9/f7xv+MhHgK9/PfCfmbnRHwA7AOzK+ezPAbwr8XoPgA0Z1/Gt\nt976/M+9997LPti5kxlg/uY3vf5tDs6dY96wgfk1rwmXwcx8++1Gl09/OlzGzIyRQRSny2c/a+T8\n7d/GyVm/3siZmQmX8fnPGxl//udxumzebOScPh0uY3jYyPj934/T5aUvjdfl8ceNjA98IE6Xt7/d\nyDl2LFzGvn3MfX3M//pfx+nya79mdHn00XAZBw8aGW94Q5wuv/d7Rs5DD4XLOHvWyHj5y+N0+Q//\nwcj5538Ol/G1r93LwK28eLHxlYYC3P112zOMOwH8IgAQ0Y0Axpg5M8bcuXPn8z9DQ0NeX/Lww8Da\ntSbiC8WjjwLLl5tI+syZcDn/9b8C69YBu3L3jZXj8GFgxQpgyRLgxIlwOV/5CvDCF0ZEIzBZyvg4\nsGxZXHbwve8B27YB994bLuPMGdMfW7cCTz4ZLuehh4A1a4AHHgiXAZiI8/LLgW9/O1zG8DCweTPw\n4INxujzyiGnTd78bp8sLXgB8+ctxujz4oLHHGBvYvRtYvx7YsydOl7vvNm2KsYHdu4HVq814c0R9\n5LvfNXP3q18Nl7F+/RC2bt2JmZmd+NCHdnr/f9Pbaj8J4NsAXkBEe4nol4nofUT0PgBg5i8BeIqI\nngTwMQD/Kk/WkSPheuzeDbznPcD3vx8u4+mngauuArZsMel9KB58EPjFX4xL6/fsAa67Drj22jij\ne+op4KabzO9QPPQQcPXVwEteEqfL448DN98c5wCeew7YuBG4/vo4wvjBD4Bf+qU4J33smCltvepV\ncf37yCPALbeY/gnFmTOm9PILvxBvA697nennmZkwGTMzJoB797tNEBaKPXuAt7zF9PPYWLicp54C\n3vCGuDH6/veBN70JIDIl51A88US8Dfzwh8CLX2z81LPP+v9/07ukbmHmzcw8wMzbmPmvmPljzPyx\nxDXvZ+bLmfnHmDk3pgs1mKkpM9Ff//q4mvRTTwGXXALs2AE880y4nH37jJOOMZbHHjOZwdatwP79\n4XKefRZ4zWvijOXZZ4FLLzWkETPRn3gCeO1r48h4716TpVx+eRxhPPGEcQD79oXXk3/0I+Cyy4CL\nLw4zXIt9+4CXvtRkcSdPhsl48kkzb6+8Mq5/n34auOYak92GZpOHDwMXXWQCjNig6eqr42zgzBlD\nOK98ZbwNXHZZnA3MzJg507QNtL0k5YzQAT182KTiO3bEOdenn44njNOngYkJEwHHkNdzzxlHtGVL\neJtmZszkevWrze9Qx3jggCmZbN4cF1099RRw441m4W9yMkzGc88B27fHE8ahQ8YRrV1r5k8ILGFs\n3270CoV1AJdeGm4DBw6Y9mzebP4OxdNPm/kfYwP795t5u317nD0mbSDUlvbuNf9/xRVxhHHwYLwN\nHDxoiPSFL4wnjBgb6BnCGB0N+79Dh4ANG4BNm4yxhNYYrTPasSM8Yjx40OixerUhj9C1kJER06aY\n6GpkxKzJrFxpCDVmom/aZEpBoTLOnDH9sXq1aVOoA3juOeNcN22KI6+DB017Nm40f4fAOulYwti3\nz8jZsEHOBkJhnfTFF8cRxubNRp+YNS9rAzFB07PPmrZs3x7npA8cMH27YUP4vBsZMfPNtie05Gdt\nIJS8eoYwQiM9OxDLlwMLFoQvEo+OmgkRM9FtdEVkFuxC23T48KyxhDrXQ4fMJAfMInzoGpGNrmKM\nZXTU6EBknGOo8e7bZ4xl3brwvp2YMEHFsmWmf0IJw7YphgCnpsy4bNwYN1+sDcRmGNYGNm4M18Xa\nQCxhWBuI6V9rAytXmsw2dHu5tYGYoOnwYTPGg4Om5Bfav0kbCAkwLnjCsNEVYCZqqMEcOWJKFGvX\nAkePhsmw5RsgzmBGRszkinEAR46YSQWETy5gNrrauDGurm11iYmkrZOOaY/NmIhkCGPNmvD5Mjpq\nsq4FC8x4h/avtQE7RtPT/jLOnTOlwosuMjYQGmBYG1ixwsg8fTpMjs0wYm1g7Vqgry9unJI2EEMY\nEjZg7To0aLrgCcNGV0Cc0VkHEGMsSScdQxg2ulq7NvyGrqNHjZEAcW2y/RtjLKOjZmyAOMM9dsz8\n//r1ceUbO18kCGPVKrNYHeKkbXsAmQxjYMA46pA5Y50rkYwNEIXbwKlTJgtcutT0T4wNrF1r/l67\nNmzOTE8bOevXt8MGrF2H2oASxsjsQKxeDRwvOtkqB1NTZpfKypVxxnLsmNEBCDcW5tk2rV4dbizW\nAQDhxgLMtsmSccgaUTK6Wr063lhsiSHkzv7R0dl+iXFGVk5/v4nKQ+ZdktRjCSPWBtLzRcIGQgM4\nm10QydlAaFZ64oQpd/f3xwWkaRsIadP0tNFn1SotSQUPRDJKW7UqzFjsJO/rkzOW0NLWqVPm99Kl\npj3HjoU56aQzCl3DYDbfv2oVsGgRsHDhrH4+sPVbID5iXLMmLgo+fty0BwifL8BshgGER4xp5xpK\nGGkbiMkwgLiy7PHjs20KnXdJUo8hDIks285/wPwOvSckbQMh/Ts2ZjLI/n4tSQUPxNjY7ICGTq6k\nsVgZIbsYkg4gdHKdPGkiaABYvNjUt0O2oUpkGKdPGxIdHDSvQ9uU7JdQY7HklSTBkDYl50uocwXk\nCSM0MwDm20BshrFmjUzQtHJlvA00nWEkCTAmwEiPdUibkgS4cmXYIn7PEMbERJiTPn58dnKFDuiR\nI7MDsXChSUFjMhUg3FhOnDBRhIXE5Ap1aMloHAjv35MnTdkGCG/P+LjJcgYGzOuVK8N2xElkGNPT\nRh/bplAHe/Ro/HwB5ttAaNCUjMZjS5hWl1gbkMowJGxgxQqTYYesVyVtIFSXZHv6+sKy/Z4hjGXL\nwu50HRuLN5ZkhGblhEz05OSKia7sxALiCMMa7kUXhTnXpPEDpk2hhGEdgISxAKZNsWMUShjj46Zk\n2N9vXoeWcJIZU2h7ZmZkCDmdeYWSV7KEI2EDK1aYDHtqyl+OtA309ZlgMkROFTZA6XPAHdAzhBE6\nuSTS8RMn5jrpFSviJ5eEcwXCHUAyrQ81FskMIzZiTJN6kxnGiROzfRuji4RznZgwh1QuWGBeS9jA\n4KA5GcC35JHOvGIIw84XIplx6kUbCEHPEEZMxBibYaTLQDHRSGyUliYMiTZJtMfqEtumFSvCMsk0\nqUtkGKHOdWxMJsBIOoClS42D9o2kk/MfkJkvRGFzJrkoC8gFTSFtmp42mcmyZeZ122xgfNxfRtoG\nQtAzhBESpU1Nma2VdlLErBuknVFo2hkbXaXJa/ny+MklEY1bObEOILQ9aSci0abBQVPS8b25LGu+\nhJBgOpIOddLpMWrKBtLlVEkbmJjwk2GP5O/reEipDKNNNhCCniIM38llIz1bywuZWICMsZw/b8hr\nyRLzWmoNI7RNSTmh0fj4+PxIL9ZYli41Z0v5LhxmZYGhJUwbkYc66XRJKjTDyCJB3zalM4wmbWB8\n3Hy/RUw0nrYBXwebbk9ogCFhAzZztLsNY4JAJYwOQgkjbSwSqV5IxDgxYZyhJS+pdDykTTYdX7rU\nvA4td0xMzGZvgFwkvXSpv1PLil4lnPSyZf66pEtSEv0CyNjAsmUyNhBCgjaqt5CygZA2SZWZs2zA\nVxdLOtY3hI5R2gZCcEETRjqiCTF+QMZYJibm6rJ4sXHcvguHEsZy8qTRxabjROEOQKJ/JUgwLSPU\nAaQjxpCIXGqTRFabYm2gyQwjbQOhtXqp+SKx5pVlA7FzV0tSAgjpxLzB9L0zuop0nMjo47tXOq1L\niAPIikRCHUAyugohDLvOZEt1QHiJIda5zszMzbyAsDalS1Ixa17JNq1YEVarj3VogEyWXVUAJ0Hq\nUhlGSJuyCMOemuwDLUklEOJc0xN04ULz4/scCqnoKjmxALnJFetcARkSDGlPOh0PlSOxtmO3oPYl\nrEYiUAnZ+cUs278WTa9hxDpXoBobGBw0Gb/v+WMSY5SeuwsWGD/lu9lCS1IJSBgLIOMAJJwr0Bxh\nZE2sUAcrHV1ZXWLT+tD5IjFGp07NzVJC5supU7NHv1iEru2kS1K+fcs8v5wkUZJaunT25FkfSJVl\nkzZA1Hs2EIKeIYwQY8lyACGTIqtMEZLtVJFhhBhLVr9Y442RE1oekzCWdMQoZbghciYn55bYQtYw\nskhdggRDZJw5Y45csfdPSOnS32+OcwnZtiyxrbYKG2iSMPQ+jAQkI0bfgUhHjCHklY6urC6xtVcJ\nAgzVJSu6koiKQowl3b8SJcxQXbLmy+SkXyQtRV7pNg0Omo0WPs9wT7cHCHOuEmXZmZn5Yx0yRlXZ\ngFRmEGoD6f71hRKGgANIR4xSzqipktSpU3PbA8g4gCajq7QDkMpIJTKMkEjankcVq0veZgsfOen2\nSOkSIsfO3XS204QN2HWmttpACJQwIico83yDCU1d21KSqsoBNGksaQcwOOh/A6BkSSorIvd10lUQ\nBuDfv1nzpSkbkHSusTZw5szsRppQGYBMm5hNQGJv/gvFBU0Y6QU/wH8gzp41C4/JxcfQ1DWWvGZm\nsne9SBiLRIZhnWJs6UWiTX19s6UgV0iWpLKckU//ZpWBmsps80pSTdhAVp1eqizrawNV7XwE/Ntk\n15n6Ij3+BU0YEsYiGV1JGO7g4Nx0PMRY8pxRbP/ahVGfrYlVRoy+Tk3KAUjMGYn2ADJBU1403oQN\nVJWRhuhSVZkZkBmjEPQMYYRGNGkHIJWOS+giYSyLFpmyi8+xHhLpuD3eOp0C+xKYpAOIJUEpIpWI\nXqXKhhLzTqI9QHUlKamyrFSG0cSid9bcDUHPEEZIRJM10X0dWtZALF5sHKbPThOp6Cotwy5i1l2T\ntsaSfkiLry5S9dvJyfnk5TtnsvpFclNBE2sYeTbg64ykyKsKGwgtSVWRYTS1tVwzjBRCbvLJmxSx\nA2EPyPNxRhLRVd4uCAkS9NUlbwuflAPwGaNz5+avM4XoIhXV5zl73zUMKV1i5bSpLCtBgEB1NhC6\n2aKKsmEIeoYwQrYm5kWMscYChE2uKtZTgGZ2vWQRINCMA8gj0pCoXsJJS2zZlNAla4cfINO/TZbH\n8vp2ZkZGjiuyCLCvz8j12WwhlQV2fUmKiG4ioj1E9AQRfTDj8yEiOkFED3Z+PlwkT8IZ+WYYeQMh\nMbmkCEMia5IgwBA5Eg4ty0GH6hJruFNTxlEnt1oCcuTlM+fOnp1/h7aVE3LvQxJNZhhpXfr7TWRf\nNyEXZdkS9nhBZRhE1A/gowBuAnA1gFuI6KqMS7/OzNd3fv63IplSzqiJDKPqklQsCUos+AHNEEbR\nGNXtROwYZa3t+OoSG9VXmZEODJjfrsfzT02ZNb/Fi+e+31RZtqpttUC7bMAXTWYYNwB4kpmfYeYp\nAH8H4C0Z11HGe5mQqklLDETI5KqyJFV36SUrWpTSRZIw6u6XvGzHd75kyVm82GQNrvVxyYw0y0n7\n9E3RJokmSLCqbbUhcpQwDLYA2Jt4va/zXhIM4OVE9BARfYmIri4SKJHWSywQW13aUpJqYg1DMh2v\nsmwY6+yXLPFbxJQKMCQ2W0gFGBIkWKVztXLaYgNNBE1SaxgLyi+pDC77mR4AsI2ZJ4nojQA+D+DK\nrAt37tyJkRHg9tuBc+eGMDQ0VCg4r5bcRPSa9WAeXxllusROrjZFV/aAvOnp+fV3VxlWl1gnbRcx\nT51yOzq6KBofGYnTxcqZmHDTJc/RN2EDkpskYttk71myZbVQXcbHgfXr57/vI2dqygQCWWtevgvn\nS5YAw8PDGB4edv/HFJokjP0AtiVeb4PJMp4HM48n/r6LiP6MiFYz87G0sJ07d+Khh4C3vhUo4QoA\n1ZYpAP9Ib9Gi+MXHojaFTK4YXaqs3xIZ0jh9Ovs7XGRYXfbtm/9+iJxYJx2SYcSWgSRtIDbLlpov\nRRscXG1Aaowksuy8MQrZabVkCTA0NDeYvu2229yFoNmS1PcAXEFEO4hoAMDPAbgzeQERbSAyVU0i\nugEAZZGFhU/EKGUseamer+HmyZDYDRQ6uZIYHPSrj1dJGIBfmyTGqEgXn+hVwrkCxY6xCRvIc7AS\nNlD3fJHcKh9rA1KE0fXbapn5PID3A/gygEcBfIqZHyOi9xHR+zqXvQPALiL6AYA/BvCuIplSAxF7\nLwcgU0v2vclHOhpJwrc+XrWTlnAAEusGVo5rmySca5EubbKBJshLar7kzd3YAyKtnLoJQ2rRu8mS\nFJj5LgB3pd77WOLvPwXwp67yJAZicNB8xjx/x4aPHN/oKkuGPVHVpz6eN7n27p3/fhbsuVOLFs3/\nzFeXbdvmv++TNU1Pm7WKLF2aiBiLsgMJXSRKUhLzTuLGMquLD2FUGWBYu3ZBHaQeW5KS8lO+6Jk7\nvQEZY+nvNwtMrvvHi6KI2OjKypEiQR8ZWZNQQhcf53r6dL4uTRFGVdGrZEmqzvlStS5S29zbNF+W\nLnXP4PL6tr/fLMq7nvyshJEBX+bOq+mFONg0JKIrKydWFwljsbrEkqDv2kOeLlJrGK7tKcp2fMZI\nctFbInrN6hepckfd80VKTp6MgQGzo9E1mCzSRSKYlLABX/QUYUg4V6D+aKRMl7qddN7Eqrv2WvUY\nhWSkEtlOlaWXugMMK6eqWv2iRX4nP1dpA/bk5zrtsWob8MUFSxhVR68S0ZWvLhIRY9UZRptIXaKE\naXWR6BffbCd9jIbVpW7CkMia8nQhqr9NF5IN+KJw0ZuIvuAg4xgzvydelXg0kepVnWG0qSQllWE0\nkY7H1pIlddmwIVsX337JynZ85WTdWDYwYEjp/Pn5R8LnyckLmo7lboKfLyNvM4Xt35gbI5uygTxd\nurUkVTYdXgjgV5F9nhN33nfexVQ12pTqSWYYbSlJtSnSa1vmFatL8j6XsrvXy/rl8GE3XfKI1Eb1\np09n362fxMyM0Tsr2/FxrqdOARs3Zn/makvMZht6bObVqzZQeYYB4MPM/PWiC4jo9+LVkEETqV5R\nSaru8ljbJmhVW1BDdMmS47M1Ma89Vhef+1PynLTrtuU6nIiVU0YYk5OmH6vayeYj58wZs+bRl1Fo\nr1sXe7xI+kiPJnQpk+ODwjUMZv5UmQCXa+pCm0pSTUWvbSKMtpSk8uT4bKGuw3Bdx7qIvOpeOyvT\nJTbA8NGl6sDLR07VNu2jC1DcNz5wunGPiH4SwIcA7Ej8DzPzdfEqyMF3UqxcKSOnynKHRPTapnQ8\nueulrD5eR/3WysnaLluFLheSY2yCSKV0yVrb8ZEjWWau0k/5wvVO778B8G8BPALA40GH9cI3utq8\nOfsz14GYmTF13sHBcBlWlzVr4nSxj9uU0KVKB+BTHy9zrsePx+li5UxOAqtWhcuou/zYpui1LvJy\nses6ovEmiFTCT+Wt7fjClTBGmfnO8suaheQEddk9U0fN1NUBnDs3W2JJwz5YZ2YmW1cfXY4eLdel\n6CY3YNYBuBBGW5xRmRORuElTyjFK6RJbkurFbKcbSd2SRZntu8CVMG4jor8E8FUAtuLLzPwP8SrI\noe7SQJGMRYvcn9lQpovLcxKKJlbyOPCytFQiuio60sNHTq9GjBJR/YXkpJsoScWuy0iSepXzxReu\nhPEeAC/oXJ8sSbWOME6fdouki5yI69EgRTKSNxzVEUkXyUjKKZs4k5PA6tXV6tImB1B3xFhHSapu\nZ9SWLLBbAgzfjL/KuesLV8L4CQAvZGaXp+Q1hr4+E9mfOVPeQWWTa2ys/PtcnXQMYUg416QuLnK2\nbo2XIaXL2rVxMorWdnx16QbyaiLDkNhu3C3ZjoQufX3uGX+ZLkeOxOniC9eq1rcBXC3zldWizslV\nFrFL6eIaXUnoUnV0ZeXUFTGePWvuXM4rC9btGIv6V2JbbZscY5vKQNZBu4S8dUT1rnNGwh6lttQC\n7oTxMgA/IKLHiWhX5+dhGRVk4ZqSV52OW13qJC+pqL4OY6mLvMrqt3U6RpvtVDnvfEpSVZdl7QPA\nZhz2VUpsKijq2wULzI/LPTd1BU1tCWx94FqSuknm66pH3QNRtWNsoiRVpeEC3UuksWMkme1cdFH2\nZ653r09NmWuydtVZXVx2ChY5I9cScdHxIlYXlweAuY51HffclAWTbfJTPnAiDGZ+RubrqkfdzqhN\nJam2RPVtKknV1S+x5QUfXSYngU2bsj9LRtJFjtH2S5U72YBZeyy6xt7LlKeLRKkOmG1TzD030v0S\nI6cJwigsSRHRA2UCXK6pE3U6o26KIiQixm7slzp0cS29uDg0ifUql7JUHUTqKqeOer+rLlZOW9Yw\npGygrpLUVUS0q+SanOS4GbTNSddV7igzFteadC9GV1VngX19ppxSVnpxmS8HD5br4uoY87ZHu+ri\n2i95dyK7ypHUJe8YDVc5krvq8sqGvnLalGGUEoaDDMdnYdUDF2dkH7NYVDOtyxm5LIT2YklKyljK\navUuukhsTbT6FF1TxwK8leOSYVS9q85VTh0BhqsuRScmWBmu9phXNgTc/VTZ2k6rCKOb1i4sXI2l\nqGZap7FILoRKOemydFzCSbucA1UkZ+FCo8PUlOm/PEiWpFzKSXn3jVgZVRMp4OaM6pgvrnIkS3US\n5NUWUi9b25HyUz7oqUe0AvWnwBKRnosuZfvH62hTf79xzmfPxulSZ8RYpwMok1MnebXJBrpJlzrW\nvFzl1EXqPug5wpBY8POp91dtuK77x+soSbnKaVv0Wgd5uUb1ZYvVEmPUa4vedQcY3aJLEyWpniMM\nqYFw2VEkZSx5C2xJOXU4APuUsKISj1T/lrWnbPHRyikbpzZFr23Spe5SXZmMsnHutZJUm0jdB06E\nQURvJ6IniOgkEY13fk7KqCALCcN1PUagjqhTSo5Ev0jJkVjbcZUjudUylgS1JJUvQ8uG2TKKdEme\niB2jiw9cM4z/BODNzLyCmZd3fkqePNwMJJg7eYdqmZxectKSuhRFjBKLsj66xEbA9gE0Em2qY2dS\nm6JXibFuUxnI7p6y2XiRnLKsKVaX5MPIYuT4wJUwDjHzYzJfWS3qdIx11G+tnDocgEvqWpeTbhOR\nFj0oS1qXOksvEmXZOoIm66TrWMdzicYl5NQdNNVNGN8jok8R0S2d8tTbiehtsV9ORDcR0Z5OueuD\nOdd8pPP5Q0R0fZnMXozq69om6WosdZBXN5K6iy6xztU+D73o2I9uLEm1xR4lg6aqS5iuujSxrfYi\nAKcBvB7Av+j8/GzMFxNRP4CPwhxseDWAW4joqtQ1PwPgcma+AsCvA7i9TK5LOl5nJN1NxlKXLm0r\nSXUTkZY9zRBwL0kVtUnitFqg3qCpTcFBN2U7PnA9fPCXZL5uDm4A8KS9OZCI/g7AWwAkS19vBvDx\njg73E9FKItrAzLkPLe3V6LWOBdU6o6tu0qUuh+byJDbXfjl2rPiaycniIz2Stfq8u56npoyuErvq\nyg4ELLOBshMTfHTpxaCp1gyDiLYR0eeIaLTz8/dElPNcNmdsAZA8tHhf572yawq/V2pnUrc5o7ZF\nNLHpeNsivTqyneST2GJ1qaN/XbKduvp3asr0Xx65uericlBfmRy7SSLvSA8roy4bkCxJuT4P4w4A\nfwPgnZ3Xv9B573UR3+36uNf0dMz8v507dwIwh7cdOjQEYChXYNvKHS73YXSbY2xLdFVmLAMDs+sC\nC3KswTUN9yHNAAAgAElEQVSqHx0tvsanTXnzqq4AI6lL3iF6kllgrA24tMelzCYx7+raJBEiZ3h4\nGMPDw8X/UABXwljHzHckXv81Ef1m8Lca7AewLfF6G0wGUXTN1s5782AJY88eoKw/6ixJSUUR3bZL\nqkjOwoUmCisqd9RF6smtiXnPXpcsScVGr3VlgW3TpWxdRnK+FB0a6CKnzlKSrw0MDQ1haGjo+c9u\nu+224n9OwXXR+ygRvZuI+oloARH9DwAczvgsxPcAXEFEO4hoAMDPAbgzdc2dAH4RAIjoRgBjResX\ngFx0VRaNTE8Xn3gLaEkqD9ZJd4vRtUkXV1KvY+NHm/qlrgBDSpe6suyyE2994UoYvwxTjjoE4CCA\nmwG8N+aLmfk8gPcD+DKARwF8ipkfI6L3EdH7Otd8CcBTRPQkgI8B+FdlcttUv61rUtS11RIob5PL\nkR4ucuqs37bNMcZG0nXNOynyahNhSLSpTetM1haL/JQPfB7RGrWNNkfuXQDuSr33sdTr9/vIbFOk\nlzxiJG/ALPGU6VJ0HLgLebmUgSSc69RU8fMEXOW0LWKU2FZ7IZakpMirzpJULCG3LfOS2lILlBAG\nEX2Qmf8jEf1JxsfMzP9GThUZLF48e75K3jlEdUU0ySNG8iLuuiYo0SyBFa0brF9fLEeiX1zl1HHv\ng6su3RTVSxwN4qJL2xxjmwKMMl0k/VQskfqgLMN4tPP7+5i7O4ngvsupViTr43mLmK6RdNn2Rhfm\ntrrEEIakkz51CliRcwpYE7oUyXEhr6IxmpmZfQhNmZw2lTva5KSL+lciY/LRRaJf2rK24+KnpMao\nNsJg5i90/pxk5k8nPyOid8qpIQupgSh6dKfrQFhd1qyJ00UiipCYXFKEUUda77K90UVOt0WvrplX\nXbrUUZKSCAJddVmyBBgbi5Nh5Uj4qaKbNKUzDNdF799xfK8VcKl3xpY7JEsvdZSkpHRpU79I1W8l\nykBl/VL2HHkLl7Gua91AYqztbsOiRwXUlWEk76QvkhNLpC7rki5yunEN440AfgbAFiL6CGZvolsO\noOSA3+ZQx0C4pON16iIV1dcRdVo5VS9ito3UXXaruIx12TEaixeb7CrviBHXrZYSJamFC40ORc9e\nr3NNcfFi06Y8vevOsqu2AemSVFmGcQBm/eJM57f9uRPAG+TUkEVdzN1rUX2vlaR8y4ZV6lLnfEk6\nxizYdR0X8qp63tnnSpTtqpNwrlZOt4x1nfboirI1jIcAPEREf8PMrc0o0qiDudsWvV6IjrFN2U5Z\nVF9nvwCzcyYrkvbR5WTBczUnJ4G1a93kTE4CK1fG6VJH+bEXbUCyJFWYYRDRZzp/PkBEu1I/D8up\nIYsLrSRVp9HV5aSloiuXMZIo1ZUdHOgzX2LX36ycvDZJZV4SNlBnvb9MjsuJt5K6tCnbcUXZttoP\ndH6L37RXJVwmRR13Ils5Vd8V2o0pcNtKUhMTcbok9YmN6qXalDdn6s52iuS42CLgZkdbHc7PLtLl\nzJny58iXybC6lG0Jt3Ly2uTqp1q1hsHMBzp/jgLY27njexGA65BzCGAbUDS5zp6VmxSxxuK6+Kgl\nqWwk76RvWpcyOXVmpEDxnKmT1K2cbiGvtuhy7lx9Jyb4wHVb7X0AFhHRFpizn94N4K/l1JBFN6Xj\nLouPLo5RyhmVyRkYmD1iJE9GXUaXvJM+T0adTrptUX3VukjYQJucdNtIvS6b9oErYRAzTwJ4G4A/\nY+abAbxITg1ZSERXbYporGPMq49LOBFm81lZCpw8DjxGF6ldL3U5xjZF9XU66apLoXVnO23TpS1z\n1xWuhAEiehnMg5O+6Pu/dUNqINpiLED1JOiaApfJaRMh17mt1kWXum4iBOqLXussA3VLgNGN5TFX\nuDr9/wnmzu7PMfNuIroMwL1yasiibCCkDLeuHSJlcqQmqGvqWrUuU1PmULaiZ0XXpYuUHC1JVa9L\nXTbQpmCy7E566ZKU6/HmXwfwdSJaTkTLmPlHAFp3Uq3F0qXA/pwledeBGBgwTivvOPC6M4yqHYCv\nLlVGjC7HtbvqEuvQXLdaAuWLu3XWpHutLGtl5D0qoM6MNHnESMw9N0uXAgcOZH/mKqPsTvpGMgwi\nupaIHgSwG8CjRPR9ImrtGoaEc5Wq1ZcZi8uWQqDYAUhEVz61zqrrwFLZjkS/nDtnnvWd97zvKnSJ\nvT+lTJduJAw7BmfPNq9L8lEBTeti5eTNmabWMP4fAP8zM29n5u0AfqvzXitRx7pBnem4lBzNdrpb\nF99SXS+VpKTkdEtZtu5+cYUrYSxh5ufXLJh5GIBgZUwWEvu+gXoi6TY56V40lrY4NMk1L5dSnURG\nWmRHPqW6OgK4OktSZXKa8A0Sdu0CV8J4moj+PRHtIKJLiOjDAJ6SU0MWbZoUbTIWCScCyBlLW0hd\nar5ItKloEbNuUi86Dtw+itelVFdX0NSW4MBn3lVtA02VpN4LYD2AfwDw9wDWAfhlOTVk0baIsRsc\nYy/qIlXvd11nkqrVL1yYXav3JXWJdby8Wn3dAUaRnOlps9a0aFG5jDrWzrox23FF2fMwBgH8BoDL\nATwMs47R+lNrJaPXC6lW3xZjqbtfklF9utxz6hSwbJm7LocPZ3/mGmAAs04tfWRM3RkpMNu/ad3r\ntiMgv01Su+oknPT58+bHZZ2pF0tSHwfwEgC7ALwRwB/JfXV1aBNzS5akqkzH6+6XhQuNgZ87Fy7D\nRReXfunvNwaeF9X7OPoqs0kpJy0x1nUHGFZOVpvqDjCK5PiSV7eVpMoqkFcx87UAQET/BcB35b66\nOpR14JYtcXKmpkxdt64by8rkuBzpUZcuIXLS/ShFXiElnHRU70MYVWdwUk7aN2vKkiMZYMSW/NpU\nHpO0ozVr4uScP298lUupzhVlGcZ5+wczny+6sE0oivQkHICdFBI3lknUb8+cqbckVXUq3aaIcWIi\n3rkCMk7NdwE+T5eJCT8baEtUn9cmKVKX9A2xutSd7biiLMO4jojGE68HE6+ZmVfIqSKHooHwcQB5\nE7QpYxkdzZeRdcdpGslTb7Nq9T5tOnRo/vszM4a8yo5rt8gj9jZF9ZK61F16KQqaJGygLf0iReoS\n5CW1ziSRBUqvXwDlj2gteWpEO2FLHFnHekg4AB8Z9vuzdJEgHh9dkqfepr93chLYtClOF0sWLuRV\nJKdtUb3EGkbdDqDqSLpNZaCmSD2WMAYGZhfJ09uTfbPA2DFyRWtPnI1FXjRSt7EUyam7ZFKmS53p\neJEcCcJwfThVmRwpZ1S3A5AsSbWFMPLa5DtfYjePWDmxNmCPH6rKBnzt0QU9TRhVTa5Tp4Dly911\nKdsO6CojlgCBah2AFGFIOOnTp9uT7dhnjdRdH68y25mYcLeBqktSkmsYddtAlesyVZSkGiEMIlpN\nRPcQ0eNE9BUiWplz3TNE9DARPUhE3/H5jioHQjKqb0uGIUGkIcZSVU1aKguUIq9Fi8ofC1wmRyIa\nlyIvibVAQOaIfx9diu6k9yVBCRuocl3GdfekD5rKMH4bwD3MfCWAf+q8zgIDGGLm65n5Bp8vqDIa\naRNhSGUYEhPUx/iL5PiuG0xMzH/fx/itnKrWMHznS5EuIfXxJM6cMWtpLkd6AMV25NqmgYHZgxPT\nqHuX1MKFhriz7v9pogwkRch5887HBlzQFGG8GeamQHR+//cF1wZtCqvSMYYQRl404nO8ed6k8CGM\nKqP6Joxl2TIZJ11lhiGpi6ucvPp4EwFG3qMC7AGGPvdhSBFy223AZ5yWLcsPmnz6xQVNEcYGZh7p\n/D0CYEPOdQzgq0T0PSL6NZ8vkNqyWaUzkko7fXQpmlyxNekmSlLLlwPj4/Pfb8JJL15sItf0YX1S\nTnp8PD5rCtFFygbSckKyHYksO8sGmGU2AzRlAxJZtgsch8ofRHQPgI0ZH/2vyRfMzESU84BBvIKZ\nDxLROgD3ENEeZr4v68KdO3c+//fQ0BCWLBmaN6AzM34DWjQQEs5ofFwm2/E1llgHW6RLm6IrH2OR\n6F97WN/k5Ny+DJkvIyPz35dw0r4yli/PPh9LwgZ8CbBojFxPbwCy58yZM6Z0JkFesTbgu86UZ9NZ\n/mV4eBjDw8PuCqZQGWEw8+vyPiOiESLayMyHiGgTgIwpCTDzwc7vUSL6HIAbAJQSBgDcfvv8gbA7\nZ1wXHyWcKyCz06QoY/J1ALEkKBUBS5WkstrjQ8ZFuoRmB7GEUVX/hkTjeTYQmzVJlZJ825SVlUrN\nl/FxYNWqODlnzxoflfVo6Cz4lKSGhoYwNDT0/OvbbrvNXVk0V5K6E8B7On+/B8Dn0xcQ0RIiWt75\neymA18McguiErMkVMrGqyjDOnzcTQ+IcnVhjAZohDIl99facnPTBgU2sYeTJaaKECciUpKq0gSYC\nDCCbBJsi9awxCskCJQJbFzRFGP8HgNcR0eMAXtt5DSLaTERf7FyzEcB9RPQDAPcD+Edm/orrF+QZ\nbt3ONU8Xmxm4nvOSVx8PMZa0A/Alr8FBk8Jn6eLrACQixqw2NbGGkSdHIhoHwhysREmqKhtoKgvM\nIsEms8DYwFaqLOuCykpSRWDmYwB+OuP9AwDe1Pn7KQA/HvodRU7aFVUbi89g5tXHT50Cduxwl7N8\nOXDixNz37AR1Ja++PkNg6UMPJSJG33UmYNZgkqd7NrGGAeRntk1kGFWXpCSI1Dcal1iAl8gw8rJj\nCRuQCAKtLr2SYVQOiZKUnVjpm3xCnEhsdAVUN7lCUteqSgy+60yAXE06y+iaKElJlOry5DSVZefZ\ngASpS2UYUgFG3SUpKyO2+uCCniUMCcNduND8nDkTJ0ciusqTI7GGEUoYaaPzddJZpB6iiwQJrlgx\nv198nntiUeUahlS5o+5NElaXqkpSEmW2kPbkRfV1l6T6+mTKj07fJSuuPZBwrkD25AopMaQH03di\nATIOVirDyHKwTaTjgBxhnDw59z3fdSZAJqrP6pdz58yd0j4Pw6kqwADktvj6zpfTp+Pvc5GYL0uX\nGl3Sd69LkGBTNuCCniUMiZIUIFPvzJIhVZJqKsPIcrASab0UqfvqktceX10kIr0i5+pLXllyYtcw\n7EO7fM4pyrMBnzGyx/PHZvwSNtDXl+2km9j6DOSXZXvlaJDKIcXcEpNLwqEB7VrDyHKwIVsKqypJ\n+RJyXoYRQhixmy0kghQpXQYHZx/1aWE3JbieBAzIlkKTbfK9yQ1onw1ohtECSOxBB2S24Emxv9S9\nJVVlGL1akvKBRJskyn2AjA0Qze/fkH6pygbsIYI+60xV2YA9XqTurc+AEkY0JBaTgGozjDbtkvI1\n3IsuyiaMurd9AnJjNDExd0dciC4SEePgoInokyfNhs6X2JIUML9/pTJ1ibJsCHlVlWGcPu13NhbQ\nLhtwQU8ThlT0mhwIe/NcbEQTEl21bQ0jfT9H29JxH10WLDBjmnZGUmsYIVF9rJOWKEkBMrpUVZaV\nIEAgjLzSQVOIjKpswAYcrk+cdEXPEkZVJamQnTNSEzSdjofWb9tUkjp1Kj6ql1jDAOa3qak1ryxd\nmipJAfNtQDLDiC1JSWYYvrqkg6ZQAqyiJGVl+PgpF/QsYUgNRFXpuMTk8j0eGjARx/T03AfISDi0\nkIimv9/sekk+J6FNTvrECfNejIxQXSTmXd52bgldmorq022S0qUpUs/btt9E2dAFFxRhnDhh0kgf\nSKbjyUhaYnKFGFzWIqZUNO677dPKSUdpsdFVSOZldUm26eRJ//kitXieNdYS0WtIm6RsQCKql3CM\nVWXZIWOUVdqVzDCk0bOEkbXTJCRilJigAwNmC2LyRNXQmmlycoUQICBTYpAwFkCmDpx2APYZ2j6Z\nF5BNGL7zJWszQFMLzXk24Dtn2rToLWEDEjsfARkbsPMlHUzGEkYV92AAPUwYdv940kmHRFcSztXK\nSRtd7M6kUMKQiBizjCWkX9IRloQDiNFFoiSVbI/vScBJObEZxkUXAWNjc98LDZqqWMNoygbsOkjy\njvGmCGPhwvmbLU6cAFau9JOjJalIEGVHI7EZRihzS8i5EDKMEGORSserKEnZMfK5yQ0wfRm7AJ+e\nL+fOmbUmn5OArS7JuXvyZFgpKXaTBCBjA/39syc/WzQdNMWSoJakBLByZfzkSkfjY2N+T9TK00Vi\nK2sMYaQdQEjpJalLKGFIZRgS0ZVEv6SNf2zMnwCz5MSUb2y5wxKg7zqThA0MDprNFjbjDzkbC5AL\nmtJtktjgEBM0Jds0Nhbvp0LJqww9TRjpgQgtSSUH4vjxMAewcuXc8kCTGUY6Gjl+3N8B5C16+yKd\nYYQYSzp6DZEBVLNLKnSM0vMlhAQXLzbkYM9eCpn/wHwbCCFBorlBU+gmiSqy7PPnTbYRa48xQVPs\nnJGad2W4YAjDpuO+teR0+SY0Yly50jhmiyYXvdNRfUibpKIriQwj7VxDCNDqEluSWrzYOCC7bTl0\nvqxaNbdNoZltcs6EECBQjQ1IReMSNmD7xbdsKBk0WV2mpkwW5usbpGygDBcMYdiJFRLRpA031gFM\nTZmor6lF79Wr55JXyORatsxEZfZ455C6NiCzhrFkiVnAtJG0JGH4Otj02llMhpEco2PHZAgjRJcq\nbEBqjCRsIJSMs+ZLrA1Y8vL1UxI27YILhjBC0/H0QMREV0ljWbky/p6FUGNZtco4IMA4WWb/IwT6\n+sxuExt5Hjtm+soXWdmOb5uITJukHUBoaSs51jHONe0AQvo3WQZqkw2EzhepDCNpA1LrTBI2ENMe\nJYxIZDlpX6xebSaCXThs0ljSC4cS0ZWdWCFHCCQNJrRNtn8B067paf+yITDXAUhEr8zA0aNznxPu\nijVrzP8CYRkTIFdiSGYHsTZgIWUDoZtHkv0iZQNShBEyX5L9G5ORnjgxu1VYCSMASdY9cgRYu9Zf\nxsCAib7top9E/TY0WiSa64wknHRoe4D5hBHqXI8cmatLCHklxzqGMGykNzlpsqgQ8lqzZm7/xkaM\nzOFtkrABScKwuoTO3eT8t3JinXRoRmo3j1gnHdOmtA34or/f6GPnb2ibytDThLF27exAHD0aZiyA\njINN1m9DJxYwt01HjgDr1oXpEutcAeNgbZuOHg1r09q1sw4g1KEBMm1KygjNLgDTD7ZNErqMjxvi\nWrjQX07aBmKdK3McCSazndAS2/j47AOdQudMer6EOuklS2aDyVC7bpMNlOGCIYwjR+IcgDWYUMeY\nLJmEpuPA/DbFRoyhxgIYsorNdpIR4+hoGAEC2WU2X6xbZ3QA4ggj2abDh4H16/1lpMupofMl2abQ\n+bJsmSkXnj1rnOPAgP/9E8B8GwiZL319s/N3aspkuLFltlgbiCXkNtlAGXqaMNIOIDbDYA53AOvX\nAyMj5m/JDCOWMEZHw9oDmP87fNj8HVOSSkZXocYi4Yys8c/MhAcGwPz+DWnT8uVma+6ZM7LzJWSM\niGad0eHDwIYNYbqkbSA2aLIy+vv9ZSTHKNSmgVkbsGMVsktKijCSNhAzf4vQ04Sxdu3c6Cq2xHDi\nhImuQuraGzbMGkvMYFpjmZ4OT+uTNdNDh+IcgCWMmMxrbMy0J8ZYkv176BCwcaO/jIEBs/NrbKz5\nDINotk2jo+HBjmTQdPSo0Sd0vkjbQEyAkbSBkZGw+QLM2kDM5pGkLhI2cOqUuRco5J6bMvQ0YaTT\nxVBjsQMRE4ls2GD+nxk4eBDYtClMjnVGto7seyIrYAx1YsJERDEOwBrL+fOmNBASMS5YYKKysbE4\nY9m0yRAFcxwJWgcbSxh23sXMmY0bTVsOHAA2bw6TIVWWlbIBSxgSNhBT79+0yegAyARNhw+H6yKV\nZVsbOHTI/C398CSgxwnDDsTMjBnQ0IHYsgXYvz/OuS5ebH7GxowD2LIlTM7ataYtMcbf1zdrMBKE\ncfCg6dsQ8gJmHWNMJL1xo9HjxAmzOOx7nLjF+vVGjxgnvWULsG/fbGkrpk2HDsU512SG0bQNJAlj\n//54GxgdDbeBTZvMGDPL2MD+/cC2bWEy1q0zpaTz5+OCJmsDoRm2C3qaMAYGTFo2Ogo8+yywfXuY\nHGssMdEVMGsw+/eHO6Pt24G9e81PqMEB5vsPHJBJx/ftCzcWwLTpuefiDNcSoI2uQmHr7M89F96m\nbdtMnxw/bhaMQ3Y3AXMJI3S+bNxoxnly0ixYh87frVvjbcCS8fS0aVesDezbZ/QKwfLlJsA5cULO\nBkJ1WbjQyDlwwPRLaP8mbaCnCIOIbiai3UQ0TUQvLrjuJiLaQ0RPENEHQ77r0kuBJ5+MiwBsxHjg\nQJwzshFATHS1YwfwzDPA008Dl1wSrosljJh0fOtW41hjjAUwDuDZZ4Ef/ciMVwikjMX27969cfPl\n0CHgqaeAiy8O18XOl5h5t2GDIYvdu017fM9LsrA2sH9/uC42gHvsMfM7ZKcVYPpUygaee85E96FZ\noKQNPPOM+Qltk7WBmIy0DE1lGLsAvBXAN/IuIKJ+AB8FcBOAqwHcQkRX+X7RpZcC3/62qdv7Hn9h\nYaOrH/4QeMELwmQAwJVXArt2mXp/aNopSRjPPRcn59JLZ2XEGMvFF8cTxrp1pl8feywu27n8cuCJ\nJ0y7QjPSgQEz377xDeCKK8J12bHDkM6+feHROJHp03vvjSMvKRu44grga1+Ly46TNrBjR7iczZvN\nGG3fHl5OvfJK4PHHZQjj/vsNkYYeS27J65ln4mygCI0QBjPvYebHSy67AcCTzPwMM08B+DsAb/H9\nLglj2b7dRHkPPAC88IXhcq66CvjsZ4HLLguP9FatMrXxBx+MI4yrrwY+9zkTxYZO0EWLzCT98pfj\nDHf7dtO3zOE16f5+06ZPfAL4sR8L1+Xyyw3pHDgQT4J33RVHGNdea8Z5927gRS8KlyNhA5ddZvpl\nz544G7j6auAzn4kjHUsYMdE4AFxzDfDJT5rfobj8ckPqjz0Wp8v27cA//ZPp5xhdRkaA++4Drrsu\nXE4R2ryGsQXA3sTrfZ33vHDNNcZwf+InwhUZHDRG8s//HG8s990HvOxl4TKIgB//cTO5rr02XM5L\nXwp861txxgIYw//a14BXvSpcxo03Al/5ipEVs7PjxS82YxRDGFdeOZsZ+D6ZLolXvtKM0ZVXhsu4\n5hrgkUdMvT30xjLAzLu77wZ+8ifDZbzoRSa7GBkJz3YAEzR985txNrBli1kgfuaZOEK+4QZTfYix\ngSVLTFbw7W+beRyKl73MBF4x/qW/3/Tvd75jfEQVqIwwiOgeItqV8fOzjiK4/JJyvPWt5vfb3x4n\n54YbjPGHlikA4BWvMGWxoaE4Xd79bkMWMVHEtdea0sl73xuny2/8hvkdM0GtU/35n4/T5U1vMnrc\ncEO4jEsuMfXsGBkA8PrXm99ve1u4jCVLgJtvBt71rjhdfuVXzG9rCyHo7weuv97MlxhStzq89rVx\nuvz8zxtZMU+Ve/WrzYLzO94RLgMAfvVXzRb3GF3sfLH2FIp3vhO46abq1jCIWcQvh3050b0AfouZ\nH8j47EYAO5n5ps7r3wEww8z/MeNavvXWW59/PTQ0hKGEVw55rkEaMzPGUKrY2+wLZrPjpYobc0Jw\n/nx4Ddji1KnwrbDSmJoy4xzbprNnwxd2pRH6wKIkmNsx/wEzRtPT4euS0pieDrvjPInJybis1gXD\nw8MYHh5+/vVtt90GZnYe1TYQxr9l5u9nfLYAwA8B/HcADgD4DoBbmPmxjGu5yXYoFApFN4KIvAij\nqW21byWivQBuBPBFIrqr8/5mIvoiADDzeQDvB/BlAI8C+FQWWSgUCoWiHjSaYUhBMwyFQqHwR1dk\nGAqFQqHoPihhKBQKhcIJShgKhUKhcIIShkKhUCicoIShUCgUCicoYSgUCoXCCUoYCoVCoXCCEoZC\noVAonKCEoVAoFAonKGEoFAqFwglKGAqFQqFwghKGQqFQKJyghKFQKBQKJyhhKBQKhcIJShgKhUKh\ncIIShkKhUCicoIShUCgUCicoYSgUCoXCCUoYCoVCoXCCEoZCoVAonKCEoVAoFAonKGEoFAqFwglK\nGAqFQqFwghKGQqFQKJyghKFQKBQKJyhhKBQKhcIJShgKhUKhcIIShkKhUCic0AhhENHNRLSbiKaJ\n6MUF1z1DRA8T0YNE9J06dVQoFArFXDSVYewC8FYA3yi5jgEMMfP1zHxD9Wq1E8PDw02rUBl6uW2A\ntq/b0evt80UjhMHMe5j5ccfLqVJlugC9PGl7uW2Atq/b0evt80Xb1zAYwFeJ6HtE9GtNK6NQKBQX\nMhZUJZiI7gGwMeOjDzHzFxzFvIKZDxLROgD3ENEeZr5PTkuFQqFQuIKYubkvJ7oXwG8x8wMO194K\nYIKZ/8+Mz5prhEKhUHQxmNm57F9ZhuGBTGWJaAmAfmYeJ6KlAF4P4Lasa30arFAoFIowNLWt9q1E\ntBfAjQC+SER3dd7fTERf7Fy2EcB9RPQDAPcD+Edm/koT+ioUCoWi4ZKUQqFQKLoHbd8lVQgiuomI\n9hDRE0T0wab1iQUR/RURjRDRrsR7q4noHiJ6nIi+QkQrm9QxBkS0jYju7dy0+QgR/ZvO+13fRiJa\nTET3E9EPiOhRIvqDzvtd37YkiKi/cyPtFzqve6Z9WTcK91j7VhLRZ4nosc4cfalv+7qWMIioH8BH\nAdwE4GoAtxDRVc1qFY07YNqTxG8DuIeZrwTwT53X3YopAL/JzNfAlCP/x86YdX0bmfkMgNcw848D\nuA7Aa4joleiBtqXwAQCPwmx5B3qrfVk3CvdS+/5vAF9i5qtg5uge+LaPmbvyB8DLANydeP3bAH67\nab0E2rUDwK7E6z0ANnT+3ghgT9M6Crb18wB+utfaCGAJgO8CuKaX2gZgK4CvAngNgC903uul9j0N\nYE3qvZ5oH4CLADyV8b5X+7o2wwCwBcDexOt9nfd6DRuYeaTz9wiADU0qIwUi2gHgepgNDT3RRiLq\n6yiFRM8AAAQpSURBVGzSGAFwLzPvRo+0rYP/DOB/ATCTeK+X2pd1o3CvtO8SAKNEdAcRPUBEf9HZ\nferVvm4mjAtutZ5NGND17SaiZQD+HsAHmHk8+Vk3t5GZZ9iUpLYCeDURvSb1ede2jYj+BYDDzPwg\ncrbCd3P7OngFM18P4I0w5dJXJT/s8vYtAPBiAH/GzC8GcAqp8pNL+7qZMPYD2JZ4vQ0my+g1jBDR\nRgAgok0ADjesTxSIaCEMWXyCmT/febun2sjMJwB8EcBL0DttezmANxPR0wA+CeC1RPQJ9E77wMwH\nO79HAXwOwA3onfbtA7CPmb/bef1ZGAI55NO+biaM7wG4goh2ENEAgJ8DcGfDOlWBOwG8p/P3e2Dq\n/l0JIiIAfwngUWb+48RHXd9GIlprd5gQ0SCA1wF4ED3QNgBg5g8x8zZmvgTAuwB8jZnfjR5pHxEt\nIaLlnb/tjcK70CPtY+ZDAPYS0ZWdt34awG4AX4BH+7r6PgwieiOAPwbQD+AvmfkPGlYpCkT0SQA/\nBWAtTD3xdwH8/wA+DWA7gGcAvJOZx5rSMQadXUPfAPAwZlPf3wHwHXR5G4noWgAfhwnC+mAyqD8k\notXo8ralQUQ/BXOkz5t7pX1EdAlMVgGY8s3fMPMf9Er7AICIfgzAfwEwAOBHAN4L4zud29fVhKFQ\nKBSK+tDNJSmFQqFQ1AglDIVCoVA4QQlDoVAoFE5QwlAoFAqFE5QwFAqFQuEEJQyFQqFQOEEJQ6FQ\nKBROUMJQKBIgojWd5yE8SEQHiWhf5+9xIvpoBd/310T0FBH9esE1r+w8v2BX3jUKRR3QG/cUihwQ\n0a0Axpn5/6rwO+6AOSr8H0quuxjmMcXXVqWLQlEGzTAUimIQABDRUOIpczuJ6ONE9I3OU9reRkR/\n1Hla211EtKBz3UuIaLhzXPbd9pC3vO/o/M/NRLSr8+S+r2ddo1A0BSUMhSIMl8A8SOjNAP4/mKeW\nXQfgNIA3dU7l/RMAb2fmn4B5muL/7iD33wN4feeY9J+tRHOFIhALmlZAoehCMIC7mHmaiB4B0MfM\nX+58tgvmqYlXwjxx76vmkF70AzjgIPtbAD5ORJ8GUFimUijqhhKGQhGGc4B5aBIRTSXen4GxKwKw\nm5lf7iOUmf8lEd0A4E0Avk9EL2HmY1JKKxQx0JKUQuEPl/WEHwJYR0Q3AubBUUR0dalgosuY+TvM\nfCuAUZin9ykUrYBmGApFMTjxO+tvYP5jLZmZp4joHQA+QkQXwdjafwbwaMF3AMB/IqIrYEjpq8z8\ncGwDFAop6LZahaJBdLbV/iMz/33JdTtgtt/qtlpFY9CSlELRLE4A+P2SG/deBfOo0NHatFIoMqAZ\nhkKhUCicoBmGQqFQKJyghKFQKBQKJyhhKBQKhcIJShgKhUKhcIIShkKhUCic8N8AiVMfW3BcaH8A\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f16e8e9aed0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot of the oscilator\n",
    "solver = ode('cvode', rhseqn, old_api=False)\n",
    "times = np.linspace(0,60,600)\n",
    "solution = solver.solve(times, initx)\n",
    "plt.plot(solution.values.t,[x[0] for x in solution.values.y])\n",
    "plt.xlabel('Time [s]')\n",
    "plt.ylabel('Position [m]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Simple Oscillator Example: Stepwise running\n",
    "When using the *solve* method, you solve over a period of time you decided before. In some problems you might want to solve and decide on the output when to stop. Then you use the *step* method. The same example as above using the step method can be solved as follows. \n",
    "\n",
    "You define the ode object selecting the cvode solver. You initialize the solver with the begin time and initial conditions using *init_step*. You compute solutions going forward with the *step* method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEPCAYAAABRHfM8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztfXeYZdVx5696Yk9PM8NkJsEwgAQKBNkYYdlqlAySjSys\naFlm5W9tOeiT12E/xTWD12uvw64tWcFaK7FeWZLXNlpkCUkgMwhFJOIIGBiYQUzsyal7Qk937R/n\nHfr2nRtOqBvem/p9X3/d773b9eqk+lXVqXMvMTMUCoVCoShDX9MKKBQKhaI7oIShUCgUCicoYSgU\nCoXCCUoYCoVCoXCCEoZCoVAonKCEoVAoFAonNEoYRPQpIhomog05nw8R0SEieqDz84G6dVQoFAqF\nwfSGv//TAP4WwP8uuOZuZr6+Jn0UCoVCkYNGIwxmvgfAgZLLqA5dFAqFQlGMtu9hMICrieghIvoK\nEV3StEIKhUJxpqLplFQZ7gewiplHieg6AF8EcFHDOikUCsUZiVYTBjMfSfx9OxF9lIgWMPP+5HVE\npDfEUigUigAws3Pav9UpKSJaSkTU+ftKAJQmCwtm7tmfm266qXEdtG3aPm1f7/34otEIg4g+B+Cl\nABYR0VYANwGYAQDM/HEArwfwW0R0CsAogDc3patCoVCc6WiUMJj5LSWffwTAR2pSR6FQKBQFaHVK\nSmEwNDTUtAqVoZfbBmj7uh293j5fUEgeq20gIu6FdigUCkWdICJwr2x6KxQKhaI9UMJQKBQKhROU\nMBQKhULhBCUMhUKhUDhBCUOhUCgUTlDCUCgUCoUTlDAUCoVC4QQlDIVCoVA4QQlDoVAoFE5QwlAo\nFAqFE5QwFAqFQuEEJQyFQqFQOEEJQ6FQKBROUMJQKBQKhROUMBQKhULhBCUMhUKhUDhBCUOhUCgU\nTlDCUCgUCoUTlDAUCoVC4QQlDIVCoVA4QQlDoVAoFE5QwlAoFAqFE5QwFAqFQuEEJQyFQqFQOEEJ\nQ6FQKBROUMJQKBQKhRMaJQwi+hQRDRPRhoJrPkREm4joISK6vE79FAqFQjGJpiOMTwO4Nu9DIno1\ngAuY+UIAvwHgY3UpplAoFIqpaJQwmPkeAAcKLrkewC2da78PYD4RLa1DN4VCoVBMRdMRRhlWANia\neL0NwMqGdFEoFIozGm0nDACg1GuW/oIjR4A//mPg5Mk4OX/6p8Dtt8fJ2LMH+Lu/A0ZG4uR85jPA\nJz8ZJ+PgQdOmsbE4Od/+NvCpT8XJOHEC+JVfAY4di5Pzox8Bn/tcnAxm4Hd+BxgejpPz7W8DX/lK\nnAwAeP/7gccfj5Px3e8Cn/hEvC7/+I/A/ffHydiwAXjXu+J1+cIXgEceiZOxaxfwtrfF6/KxjwGb\nNsXJGB0F3vve+PX4wQ8C69cH/jMzN/oD4DwAG3I++zsAb0683ghgacZ1fNNNNz37c9ddd7EP/uiP\nmAHm73zH69+m4NQpI2PZsnAZzMwf+pCRc+ut4TLGx40MIE6Xz3/eyPj7v4+TI6HLnXcaGR/8YJyc\n5z7XyBkfD5dx771Gxl/+ZZwu11xj5Bw5Ei7jwQeNjN///Thdfv7njZytW8NlbNpkZNx4Y5wuv/mb\nRs7mzeEy9u83Mm64IU6X97/fyHnssXAZJ08aGZddFqfLn/xJvJ3693+/i4GbePZsYysNBbjb67ZH\nGLcB+FUAIKKrABxk5ky/bt26dc/+DA0NeX3Jhg3AkiXAd74TrujDDwNr1wL798d5AA8+CCxYYH6H\nYts2YPFiYN48E7GE4rvfNf2yIbeGrRwjI8Ds2aZNzzwTLueOO4DzzgO+//1wGRMTwNatwNlnx3l7\nd9wBXHwx8PWvh8tgBh59FFi40PwOxb33Aueea6KVGDzyCLBiRdxYb9hg5t1DD8Xpcs89wEUXAd/6\nVrgMu44kdLnwQuCrXw2X8dBDpl+eecaMeygeeABYtiwiOgCwePEQVq5ch4mJdXjf+9Z5/3/TZbWf\nA/AdAM8hoq1E9GtE9A4iegcAMPNXAGwmoicBfBzAb+fJ2rcvXI9HHwVuvDEulN6yBXj+84Fzzokz\njA88ALz97eZ3KDZvBp7zHENgW7aEy7n/fuBXfzXOoG3aBFxwAfATP2HSQaF44gngda8DnnwyXMbm\nzcZAv/zlwH33hcvZsgV41auAp58Ol7FjhyGw666L65cnnwRuuMH0TyhGR40+r31tnJxHHwXe+lbz\nO9RpYjb9e911ZrxC8eCDwBvfCOzcaVLOodi8GXjNa+J0eeAB4NWvNm2LceA2bjT9G0PqmzYBV1xh\nnIMf/9j//5uuknoLMy9n5pnMvIqZP8XMH2fmjyeueSczX8DMlzJzrkkPzeGOjZmF/7KXGc88FD/+\nsfH01q4FnnoqXM5TTwHXXx8v4/zzjT4x5LVtG/BzPwc89li4jMcfN97iqlVx/bt1K3DNNXGEsXEj\ncMklpm9i+mXLFmBoyMiYmAiT8cwzwJo1wPOeF9e/mzYBV11l9nZC9722bwdWrjRORsxeyGOPAZdd\nBixaZHL/Idi7F5g1C7j00jgj/eSTJgpcuzZ8zhw/bgz81VeHGVeLbduA1avN3At1vsbHzbq+5hoz\nXqF45hmzFi+4IKxf2p6SckaoJz08bLzO88+PGwhLGOefHz7RR0fNxvvFFxvPKBSbN5uFsnp1+ERn\nNl7nlVeaPhofD5OzZYvpk5Ur4wjjmWeAF73IEPyBokLsAuzcCSxfbqLAHTvCdXn6aTNGZ58dPk47\ndhhdli+PG+unnjIpk1WrDKmGYOtWMz5r18YZ6W3bjB7Ll4f3r8Q6AgxhnXOO+Qklr6efNmso1sHY\nvt149OeeG74Gdu0C5s83pB67jlavNoQR4pT2DGHs3h32f8PDwNKlZpJv3x6eY7QTPWbh7tplcpQL\nF5ow+sSJMDl24Z57bjhh7N8P9PcDZ51lJurevWFydu82bYohjBMnTMXWsmVxBmB42MiIMWh2H+Tc\nc82eSmhaypLXsmXh7QEmDeOqVeFGbds2Mz5Ll4avI2ByLa1YEe582XW0enWckd65c3K+hBKyXUcx\njhcwSRhLloRX1u3ePbVvQ+3U1q2mPaGOihJGZ5LPnQvMnGkMUwisAVi8ONy42kne12d0Cp1ce/YY\nPVasCDeMdpIDcbrs3m0WSgxhbN9uFn5fn5EVOtaWkGO8+kOHDJHOnm36JTQnbSOMGMKYmDDEvnBh\nXMrPEkZM3wIyhCGxjqyc2AjDrqNFi0y6L7Sk266lmHU0PGzGp78fGBgI7xs71osXh81dJYzOJAcm\no4wQ7N1rJlboQACTkxwwhiTUqCUnemgxgDVoQLxntGSJzMIF4vrXEkZMSsqOM2B+hy7cHTuMHjFG\nZP9+YHAQmDEjrl9s5GVlhHivJ0+aqHjBgrh1tG+f6deBAUOIo6P+Mpgn15LEOiIypBy7lmIiOLuO\ngLioybZp8eIwXZQwEoQRY2D37TOTSsKgAXGTYu9eo0fMJLcyAJkIY8ECY+BCYPsWiO/fpUvjvHop\nwrDzbtEiE9WGVBXt2TNpRGLG2vZvf7/ZcD582F/G7t1mbPr6zO+YebdokTHSixaFjfWhQ4ZEBwbi\n11FyrEPaNDFhxnfhwvh1ZO3UwoXxa0kjDAHCCDVqp04Z72r+/DgjYr0rII689uwx/x+jy4EDpj+A\nOM/IhtJShLFkSThhWPKaM8cs5OPH/WUkjcjChfH9O22akRPSJmukrS6h/bt/f/xYS6wj4HTnIKR/\n0+sodIySkW0oIR86ZNLd06fLROpAeP+Ojxtn4Oyzw9fRGU8Ye/dOTtDQgdi/3wyC9a5CDdqBA0YO\nYH6HVAMdO2a81cHByUkekmJIGpFQT495MlKxRjokD5w2IqFjbY00UXj/SkUYyf4N1SVtRGIjDCDc\nqCVlLFgQXsmW7t+QeZdeR6H7kknCCB3rpOMVsy8j4RwcOGCKWKZN0wgjeIIePDg5uUIJI006Bw+G\nlaFKEIY10ERmY3bGDODoUX85lgStLocO+csYGTGFBDNnGn1C+1fC62Q24zJ/vnkdatSkCENirJP9\nEpOSkiAviXUETI0OQo1aeu6G2gYbqQPh/Zvs2/nzw8nr0KGp/SuhS0gVZs8QxtGjYYeo0kYklDDs\nxJo+3Xj3IRNDwogkJwXQ7EQ/fNh4NBYShDF/fjx5AaZ/Y52DUMJglvGCDx82t38BZPYwgPCxliIM\niYhfYh1ZOXYNSKyjwUGziX/qlL+c5FqSWEd9fWEHPXuGMAYGwjbr0hNdwkjPmxdm1CQm+qFDk0YE\nkCOMEF2yCCNEl+REnzcvnowBGU861LgeOWI2mGfMmJQT27+h42zz2tZpCiVkCccLOH3eha4jCSN9\n6NBkm0L3q5Lt6esz4xVip9JjHUsYgIn6fdEzhBG6eA8ciJ/oyYkFhBOGRCidJowYXWIjjEOHphJG\nqFFLGvsYI5IkjJixtv0rMc5AXISR9DpD58vgoMlrAzLraGDAlNn6pjzGx41xHxw0ryWcg76+cDkH\nD06O9fz5YYY+7UxKrKVQxyu9BkLQU4QR6xmFpinShjEmwpDIJafJK3SiS6ZMYnRJyolZ/Ml+kSBk\niXEGZCKMWbNMqsvXSCcjJiCuf62c0KKCw4cNWVjPVyJSB8L6d2LCpLkteZ11loxzIDHWobYuvR5D\n0FOE4TvRx8ZMaeXcueZ1TLiYHAgJL1gqwgid6EeOTJ2gEnsYg4Nhdw5NypEyIhL9O2eOmUO+ZyiS\n3niMLskxIgob6/QYxexhJNsUspakomOJsT5yxERKNvIKdXaSY2R1iXW+QtdR2rENwRlNGHaCWo9m\ncDCsokhiop88afKs/f3mdUzoKuHVHzky6V01TRjJNg0MGC/a9+mIaU9aYtGFGumk5wrIpKSAsLGW\nIow0CYaspaz2hBBGeg2E2oY0AYY6XsmxDnVsT56ctA1z54Y7XhphdBAyEGmvSIq5Q8J6a0Qsec2d\nG1bFIEFeExPmuwcGzOv+fpNf9j3oljYAMRM9aaRD2pQeoyadg7QRCTWM6f4N8eqrijBCxjrL0If0\nS5qQQ8Y6uX8BhI/RkSOTGQwgLCVlo5SkYxu7jkJxRhNG2ohIMXfI5Dp6dOrEmjPHGGjf8xzpiR7i\nGY2OGpKw4bj1pEMMQGyEYTdC04suxEhLjLVEBJce67lzw8grizBiU1KhexgS0aRUSirdvxK6hKar\n0+TV1DqycjTC6CDEi0h7elaG78noKrxOIuPh+0YZEgYtrQsQZtTSRBoy0e3i70vM1BCjZjdUk7r4\ntiddxWN1iR1rKcKQSEmFRl55aylGl1DykuhfqcpHCV2yIvUQO6URRgIDA/GEMWOGOXjnm3qRqJJK\ne0WAzOQK8TrTYbSULqGeXnqSDw76E2nWwvXVxfZLmrxix1oywogljFBdJPq3qghDYu7OmWP2EXwL\nHNJrSUKX6dNNVZzvbXZ00zuBkIHI86R9J3rakw7Zf8gijBADOzIyVU5oykQqwoj1XrM26iScg1Bd\n0v0SSsixXqeNduw+U6guWXPXVxfm7H2DWOegv98YaN9Dd1kpKd82pfvWpmVjo3WJSB0I61/d9E5A\nijAkUluhukh49aOjxhuykDBoobpIRBhZYbTEWIc4BmkjYuXEOgchc+7oUaNLOtoJKWVN9q8lY590\nx7Fjk9G5RWibJNKyEmM9MjJ1HQFhhJFuk8Q6snIk1pIvlDAyCCNkciUNSYguUl59mjBCF1xbUlJS\nukhEGFmEIWXQJIzI4GBYSiqpy8yZptjB5wCgVKQu0b/2oUux6zG9jkJ0AWRSUnlp2ZD9wPQ4+UIJ\noyVGOm8PQ4K8Yg0aEGZg0znppowIcLphlPI6Q52D5Fj39xsD7VMRl0UYofMudg1IOl6x/WtlJCOv\nkDFKryPAPxXKXK1zING/vlDCiJzo9kE89lBNjC5ZexhNkVdbUlJZnl5o/6ZTLyMjfqkXSa8zK/Xi\na6Sz0hS+uuT1r4+cqhyvEF2k9gIldDlxwqTp7E0mrQwpwvAlr2PHptqpEChhRBLG8ePmuRNJj6ZJ\nI532yCU2iEN1ySMMXyOdtW8Q2yZbaeLzzGgpXSQq4qQiDAlCloow8qJJH12kqg0lnAPJ1G7spvfx\n42a+27NVoegZwpDw9AB/5q4qTWHl+EY7aS9iYMBMfh8jLTHRmWXy41L9K0GCWbpIzbsmCSPWSEul\nMPPGuoloRyIlJaWLxB6GRDoK6CHCaKqsVjJlEpuSsl5EMtqZPt38+BhpiWhndNTokgzHAf+JLuHp\nMWcvGCldmhjrPMKQMtK+hlFi/02if9uUkqoqtWvlxM7dEJzRhJG+5QQQZkTSnsjs2f7141JGOq0L\n0IxnlFfCJ2GkfXU5edJENmnykujf0LLaKvpXcg/Dt1+kjHRsyk8qJZUXYTSVkopN+SlhpCC1iSmR\nkgqpH5dISeWFnU1M9Kz2AGGhdKwRySNSibA+pHJGag8jbUSa2sOQirLblJKqam+nv9//HnHpw7iA\nzDoKQc8QxrRp/sfl8yaoBHM3YaTzdGkilC7SRSLC8GlPkS4ShtFHl2PHzF5O8pBbiC5ZRiSEMCQI\nWSLdB1SXkrJjVHdFXBZh9PUZub7FFrGObU9EGER0LRFtJKJNRPTujM+HiOgQET3Q+flAkTwJAxAy\nEFnMLeGRh0Q7bUlJ5U3QJvYNinSJHeuQvs2KvCT6NyTaySq1lCBk6xjEGmmJslrrTPoYaYmUVF6U\nLWWnzqiUFBFNA/BhANcCuATAW4jo4oxL72bmyzs/f1IkU8KQSFUfhBBGFRvwQLuiHYk9oqaiHYmU\nSd5p2yb2U7KKJEJ1SffLjBnmx+dGnlVVJgHNpNmq1EViHYWgyQjjSgBPMvPTzDwG4PMAXptxHbkK\nlJoU3Wyk86KdtqSkmihbziN1ibBeyqBJeIy22MI1Py5VJCHhHGQdgAVkIgxfXQCZEuqqyetMK6td\nAWBr4vW2zntJMICrieghIvoKEV1SJNBncjFXmxtsKiWVR15NpKQkNpqr3sOI3Tj0PedSZZqCyLzn\n2jdS0bEEYUhFO1VGcE2mpNK6NJWSml5+SWVwWWL3A1jFzKNEdB2ALwK4KOvCdevWYXgY+OhHgRMn\nhjA0NFQoeGzMLLDYcwJF+wY+5CVRJdW2aKeqTe/kHVXJIf6sUpdp08wm9rFjbguyyOvcti1OF2Cy\nb1zuSipZDBC7jyexjgCZPSIb7cyeHS7D6rJ8eZwuY2NGnyw7FbIe169fj/Xr17v/YwpNEsZ2AKsS\nr1fBRBnPgpmPJP6+nYg+SkQLmHl/Wti6devw0EPADTcAJVwBoNrKGV85x4+ffs+ZUF1iU1L2kFtV\nhDFnjn8lW7pNM2caojh50niloboMDgI7drjrUkbILoRRZYSR1CVGRhPOQZWROuBnYG0hQDra8SWv\nPBL0aZN1RNKOkW+lldVlaGiqM33zzTe7C0GzKakfAriQiM4jopkA3gTgtuQFRLSUyHQVEV0JgLLI\nwsJnIPKMq69Bk5joUkZEIiU1NmYWSpq8Zs40no7rE8eKCCO2pBCQGeuQ/o0l5KNHZXSRKKGW6pci\n8nJdS2URkyuKIrhY8vIlDIk9IinHq+urpJj5FIB3AvgagEcBfIGZHyOidxDROzqXvR7ABiJ6EMDf\nAHhzkUypgYi9IZ3VxceIZBHG7NnmtLiPkY7VJa9fiPyNUdWEEatLiAGINWqS5BXbpqI9DIl+8Rlr\nCTIukxNrG6TGSEIXqXXkiyZTUmDm2wHcnnrv44m/PwLgI67yJAaiv39yE9MlPy41KbImud3EHB11\ne7TiyEi+d7Vnj5suRakV26b588vljI4CCxac/n5oKJ2GlHMgtXkukQaqe7+q6pSUz1hLRRhFTlMs\neUk5GL57O3l26tgxv328bi+rFYcEYUyfbjYyY1MvEkYEkFt0ErpIhPU+7cnbfLS6xI61jxEBqiUv\nX12kop22E4ZktNOWiFRiHfX1mTSx601Fe6GsVhwSzA34T/SqjIiVI6FL3eSV178+MrKeNWLhY0ik\n7rElRch5e2dNpOry5otEjb+N1mN0aYK88trje1NRCeegamfSFz1FGE0MRJV5SqtLnXX1kuQVu1lX\npotEtOPatxMTxpurMtpx1WV8PL9CTGIPQyrakVhHvo+vlTLSeSniup0DJYwKIRF2Av6eUZVpoCZS\nUnm5Ton+ldgIBWTIy1dGVqmllVNnv+SVWkrr4noYscqUVHIfL1YXifUolX5sIuKX2MMo3PQmoi85\nyNjPzDfGqxKPppi7ypSUlC5tIS+pMZIgQcm9nTqrXqo2aPYwYtatOtKwB8tmzjz9szlz3IstihwV\nO05ZBR1JFEWBUsa17krBtkUYZVVSzwXwH5F9PyfuvO9cxVQ12pSrl6hBt3Lakh6TSkm1ZbFIep2x\ncmbNmrwPVNlzl8vIa98+d10WLsz+zPZvGWEURTsS68hHjj1wF6uL1FhL7ON1G2F8gJnvLrqAiP44\nXg0ZtGkgJHVpiyfdaympOXOMFz0xkZ1qctVFol+SqZcyT7rMoD3zjJsuRUbazpk8QnHRRSI69pFT\nteNl5bg4X2NjJqWXPgBrZTQR8VdeVsvMXygT4HJNXehVwoid6HV7ekD7iDRLTl+fSV+4RIJ1eHpS\nhrFOY9QthCHheAHumQOri0S0E9svQM1ltUT0k0R0a+chRhs6Pw/Hf70sJDd3Y73gXiQviUqT5MHI\nGF3qJp46dHHt36rHCKifMCQclTNtvrQxJWXxWQB/COBHACbiv7Ya+HbgypVxcmx5Y+wm2+gosGhR\n9mcSBsDq4nIqtGyCxu6FzJhhPPuxseyN0iRcUiYxulg5Lv1bdXoMqN8YSfRvtxCpb78sWxYnpyxK\n8Zm7ebq4VnPmPWskBK6EsYeZbyu/rFnUfQ6jjk02n1A6zwDYqpe8CpK0LkUbobt3u+ni0r9lhFFH\nON4mw9hE6qVqXSQKP+rew2hT5CVlp/IOwPrClTBuJqJPArgTwMnOe8zM/xqvghyaWCx5E2vWLBN9\nxFa9zJnjVvVy6pT5ybvdt22TC2GsXp39meuiKzpYltSl7J5UkgagLV591fOubq9eas+rrKig14i0\nLQ6GL1wJ40YAz+lcn0xJtYowfG7IJTHRi2TYqpdjx7LvRJtE2cLdujX7sySKyhuByTZl3RDQVRfX\nReeqSxnq2sSMTUlJGum6ySu28qtbDGPd5FVkG5Kn12OdybYSxk8AeC6z69nPZjBtmvFqXfJ1ZQNx\n8GD595UNhF10MYThunDLqiDqXHQu/eK66PL6zlVG3hPLLOo00hMT+XtePnIkvfpuiHaaSElVuZ9C\nNOncltmGsjHau9dNF4mSWsD91iDfAXCJzFdWizo9o7KBqNtIV62LxGKR1EUq2qlrvhTtefnI6RbD\nKOnVx46RdRhOnsz+PK1LW5wmif6VKqkF3COMFwN4kIi2ALA31GVmfqGMGnKoc7Ou27z62EilF6Md\nn5RU3jNJJAya1aVNaaA26XL4cJyMpD51FFu0aaybSEldK/N11UPCM3ItV6vLk5bwItpkpF37t6y6\nqe6Fe8452Z9JeXo+/ZtXhu0qg7k4dSsx1pKFH7t2xeli5YyOxhVbtCnil9LFB06EwcxPy3xd9ejF\nlJRre3oxJVVlmsJXTtVep4ScZOqlyJM+dswY9LxSy4EBt+q80VHgrLOyP/Mt/KjDq5fYI5Kad92Y\nkircwyCi+8sEuFxTJ+rcrJPyGKuu2ALkctK9tlgkjEiyOq8IUgZNon8liTRWl/Hx4jNCdadl20Re\n3ZaSupiINpRc4/C06frQpoFoU9hZJ5FK9YsLkZaVUJf1y8CA2y24i9rkWp3n0i/btsXpYuWUpV7q\nmLuucuooBvCRI7XpXdT/dWdCaiMMBxmODyysBy4D4XqwrAwSRprZXJNnaNpWmdSWlJTrMxvq0CUp\npy5dJCKMMiKViDBc9qvqSKcCbvO36PkePrqMjADLl+d/7iKn7JYerdvD6Ka9CwuXSSF1sEwi3XHi\nhMk7T88ZCVdPpI5Nb9dNzLoNY4yRdjWMrv1bdDvwNhFGHfOljbqUjXUdh04BN9tQ9Ex7H11q28Po\nRrhMijYt3Lo94CLY24vkeVfJA0dV62LlxBp7KcKoY6yldHExRpLpm9ixbtN6rGMv0OrSFjvlg54k\njDZ5NLGLxRrvsbH6dCnaE6jbGPUSIbfJMNaVkpLQpU390ibb0ERKqucIo9e8KyunbPG6GFcXj6bs\nlip1e0ZlRq0ur75t5NVNaaC26FKnbWiTLrWnpIjol4hoExEdJqIjnR+Hs5f1Q8KgSZVJShJGHZ6R\niydSl2GcmCi/HbvEWPdiSkqqX9riSdeZBiozrq4PAHNZj7FjlNxTjJHjA9cI4y8AXM/MZzHzYOcn\n57hOs5CYoH19k2WSMXJcdXHx6uvw0lwmVl2eke2XovRYm8L6boswJIxrnbq49kvRWpJwmuwDwMru\nSSWhS9maTh6MLNOlbsLYxcyPyXxltZAwaEC7vHrXNsWmb2yFSBGkUlKxG+eAm0deljZ0kVFWDADU\nl6orK8N2ldOmXH3ZGNnnrk+UPOtTok0u6Zu6oiYpO1XWvz5wJYwfEtEXiOgtnfTULxHRDbFfTkTX\nEtHGTrrr3TnXfKjz+UNEdHmZTImFa+W0xauvyxhJEmm36OJTaVV1tOOiy4kTpgQ7rwzbVReX1MuJ\nE25Guuqx7uszpBEb8Us4XkB9EVxda8AHrjcfnAfgGIBXpd4PfoASEU0D8GEArwCwHcAPiOi2ZCRD\nRK8GcAEzX0hEPwXgYwCuKpLbtlx9L5JXHbq4eHp1psekosBu0cWWUI+O5t8Hyh4sK9tnkkiZ2L4p\nuq5N67GbnCYfuN588D/IfN0UXAngSXs4kIg+D+C1AJKpr+sB3NLR4ftENJ+IljLzcJ5QqVCvrhOq\nbUuPuRjpNkU7UumxiYn8A1J1pSlmzDB6jI3lP/DJtV/27y/XJe+Ot0k5RYRRdrDMyih7GJlLyqSs\nf23ZeV6/uciwuvRaSsrFyXCFa5XUKiK6lYj2dH7+hYhWRn73CgDJh49u67xXdk3h9/aiV1+n99om\n76qOKNCeyX0+AAAgAElEQVSmO4q84Lr6xW5iFsmp04iUOQd1jZGLHMl1FBvB1VXhZ+VIkKArXFNS\nnwbwWQBv7Lx+a+e9V0Z8t+vjXtOZ48z/W7duHQBg505geHgIwFCuwLYZRglPus4zIW3pF4loJ6lP\nXv+5GsayZzb4GLW8hzVJzpfYsZb0xvOeNeIqp0394hp5NbEe169fj/Xr1xf/QwFcCWMxM3868foz\nRPR7wd9qsB3AqsTrVTARRNE1KzvvnQZLGBs3AnffXfzF3ZYGqtNIl5X41pWScl24ZXea9SGexYvj\ndGmLJ91tukikpCT7JW8eSOpSd0rKyhkaGsLQ0NCzn918883F/5yCa5XUPiJ6GxFNI6LpRPQrABwe\nP16IHwK4kIjOI6KZAN4E4LbUNbcB+FUAIKKrABws2r8A6gv1XEst60pJSZ30rsMYuRyMlEx3SKRe\nes1IS0RwvdovdZGXy3os06Vsr9U+a6TMEXSFK2H8Gkw6aheAnQDeAODtMV/MzKcAvBPA1wA8CuAL\nzPwYEb2DiN7RueYrADYT0ZMAPg7gt8vkShhXwHzukteuutTSyimaXLY2v0iOyz2p6jIALgcjpVJS\nEtFBnbn6Oo10m7z6tuxhSM2XtvRL2bNGfOHziNZfkPnKKXJvB3B76r2Pp16/00dmmyZo8jYCeQPm\nqktR6mVszBjhogoRK6csP150e26g3KC5kFdSlzzPp840UFmbejEl1TbD2E1efZmcNqWkXGT4oJAw\niOjdzPznRPS3GR8zM79LThUZzJ5tQrCiZzZIDIRL3nX6dGPEiyompHRxmRQuhBGry4kTJpopel5G\nUk4eQUlFGG3z6tviSbepX9pEXnXpMnPmZFo77xCm1HyRKqkFyiOMRzu/78PU6iSCe5VTrXB58Hyd\nzG0XbwxhlBkA10nhEr7W1S8uXtpZJXcrK9PFpbzRRY5USsqmB4qgKanqdZGI+GN1SZZQ581ziahJ\nsqQWKH/i3pc6f44y8z8lPyOiN8qpIQvbibGEUZQG8jGMo6PAggVxusQakaQuMXIkDJrVpWiPaGQE\nWLasXIZE/lYqJVUkY2zM6OGaNsxDnZ60i6NSZxooVpfkjQPzHs/sqsu+fXG6AJP9G0sYBw7E6+IK\n103v9zq+1wpIeORlBk3Kq68rPVaXLt1IXhJyJIwrUG+uXsJRiV1HgAwhS/WvVErKpSpJIjpoVUqK\niK4D8GoAK4joQ5g8RDcIoOQZcM2hjkXnu2+Qh7rTQN1EGHV4nS5yRkaA1avr06XMMLre0iMPLveA\nsnIOHSrWpa6U1MAAsGNHnIykPmefHadL7JpO6pIHKWeytpQUgB0w+xev7fy2hHEYQOzBvcog4Y20\nyUi3yXuVTEnV0S8u3lUd+wa+aYoYObNmmRRY3obq8ePGAy46iWx12bkzThebgiu7P1ZdzoFU9qEt\nYy1lG1xRtofxEICHiOizzNzaiCINqdyrhDGqK9px1aVqwygVedWZphgYiN+vkkxT7C04EuvSL8nC\nj8HBMBlWl7J5l+epp+WMjmZX57kcgE3KyEPRXkBaTh1pIKlosi5dXFHoYxDR/+38eT8RbUj9PCyn\nhizaFOq1KQ3kQqSxT/+rMyU1e3bxIyolyatMl6QnnYU6jQhQPNZ1RoFlcqyMssKEOnTxPUNUpS6u\nclq1hwHgdzu/xQ/tVYmiRefyxDIroy2ToszQ17np3aaUVLI0McuT9jGukm3K8qTrTFMkdYmRUcdY\nt4m8XB5OZWWU9cuSJeW6FI21FHlJ72EURhjMbLeZ9gDY2jnxPQvAC5FzE8A2oGggjh93O1gmefYh\nb3L5nBNoi5FOPrMhVpeqiwrq9uq7xTDWmU4tk1OnswO0L/LKm3d1lmH7wLWs9h4As4hoBcy9n94G\n4DNyasiijoUrYdBczwlYGXk366szHC978Hzdi07CANQRwTVBXnly6kynlunSJiPdJiKV1EUyJeVK\nGMTMowBuAPBRZn4DgOfLqSGLtk3QWCMyY4Yx1HlevYRBc/VogPq8tNgIrhtTL3Wc6u/Gfqkr2uk2\nZ6dVKakkiOjFMA9O+rLv/9aNthm0WF2A6tvkWjtu5XQLIddVPeaiS5v6pS2GUVNS7dDFFa5G/z/B\nnOy+lZkfIaK1AO6SU0MW3ZKS8hnMqieXry5VhtLj4ybiybt9QxJFXnDdKam2GaO2pDvaRF5lqbrY\ncfbVJbZf+vuN3rHpale43t78bgB3E9EgEc1l5qcAtO5OtRYDA8D2nC1510ledjfJNhnpOnPsQPWG\n0cpwuYd/1akX170dq4sEqVcd7bg6O2daSspVhpRtGBgAtm6N02XGDFPAk3d/rEb2MIjoBUT0AIBH\nADxKRPcRUav3MGInRbJkM0ZOXV59XeG41aUboh2JA2o+ezt17WF0Y+qljminLl3KCj98bg0Su46s\nnNg14ArXlNT/AvD7zLyamVcD+IPOe61Er3r1saF0HRO0Tg8YqD4l5dq3QPX94noq2kUX1yi7qIRa\nIi3rKsMS9smT2Z9LjXU3rgEpG+MCV8KYw8zP7lkw83oAgoGOLCSMCNA+r74t0U7VpaxN9MuxY8Y4\ntkGXvJy09VzrStUVRdmuB2CldJGSI+k0Fclx6Ze61kATZbVbiOi/ENF5RLSGiD4AYLOcGrJo2wTt\nFl1cHxRfR0rKdZJLOAf2GeNZKYa6x8g+syHLq5cyaBLe68mTbqeii2QAfmNdRF4uD6eqQxcr50xP\nSb0dwBIA/wrgXwAsBvBrcmrIQmogyjyAtmw0150e68VwPK9/605JFcmpe4yA/H5pwqDlyTl+3BB+\n2d13ge5JSTUx1i4oex5GP4DfBHABgIdh9jFaf9faOibFmVxWK9G/kv2S98SxkL2QxYvjdDl4MPuz\nkLGePz9OhpQnnTXvJHVZtcpNjhR5VenVj48XP9EviTalZV1Rxsm3AHgRgA0ArgPwV3JfXR3axNx1\n6VJ3hNENKamRkfzH9Gbpk2cY2xJh1O3sFMmpO1KX0qWsX2J18d1nqjoT4tMmF5RlIC9m5hcAABF9\nAsAP5L66OpTlb885x11O1kDY/LJEtUq3Rjt5z2xw3Qupq198iKdKI113GkjKe80j5KNH/ci4LeTV\npjEq08XlWSNFck6dcj8A64qyCOOU/YOZTxVd2CYULRaJid4mI+Kz4WefH5FVDeRza5Ayj8alf4vu\neuubSioa69hIRdIAxBYV+MqoOtqJTWtJ6dKt0U7Vdson2nFFWYTxQiI6knjdn3jNzOzwjKv6IRl2\nVm1EfOQMD5/+vn02s8uGH5ExOMeOnd4Ho6Px6RtmdyOdPPyUPhQnlR6TGuvYKMXKkTDSEl69hGH0\nIeOqdZGcL7GEHKIL8+lGXSI95tO3rih7HsY0Zh5M/ExP/N1KsgBMqog523uVmOg+A5G8jUCWnDrz\n2lZOVSRoH0Djciq6SI7UvoGPge2GlJSv1ylRVpsXeUmRV5tSUhJt8tGl6C7UdWdCXNHaO87Gokov\nzWcgig4/SXivR49mP20uD1WG0j6TvEgXCYN28qRxGlz2mYrk+Bq0Kg+FtanEt026hMjIOhgpQcg+\nqV2g2nJun35xxRlJGHUPRJWhtE/EVLUuvhUZVUY7UrqMjLgTcp4uPvtMRXKkNpol+lcqJVV3hGEP\nG2bdYqTuCMPKqaqooPaUVFUgogVEdAcRPUFEXyei+TnXPU1EDxPRA0R0r893VFnd0YRXL2FEukUX\niU1vX8KQSgNlyXB9LHCZLk3sYVSdkvJZS91ipCWcSQnHtpdSUu8BcAczXwTgG53XWWAAQ8x8OTNf\n6fMFEgNRdeqlTUbapwJHykgXLVwfIxJr0IBqHQwpIyLh1bs+R74OXawcibMybelfX12q3K/ydWxd\n0BRhXA9zKBCd379YcG1QUZhUWB9rRIp0adO+gVR6zEeXuXOrjQKlyKsb0xT9/SaySZdQW6/TtdRS\nItqx5dzj41PfHx/3f9Jjm4x0Xr+4Ol5WTmx6d+5co3saR474rUcXNEUYS5nZFokOA1iacx0DuJOI\nfkhEv+7zBUWhdOxAtMlIh3j1VS4WH10GB+P7145zehOzTSkpn+itSI6Pke7rM4b6+PGp7x854p9O\njZ13eYUfdv67lIRbXaoi5IkJmbNIUulqn3knZadc4PTEvRAQ0R0AlmV89P7kC2ZmIsp5wCB+mpl3\nEtFiAHcQ0UZmvifrwnXr1j3799DQEObMGTptICYm/CZXHYRR50YokD/RfQyJlC5z55rvzZLjqkuy\nNDFZERVCGLFe/ezZJt0zMTHVCPoakSJjFOIcJOe71Nz11WVw8PQ5JrnntWRJnJxjx8z4ue4z5eni\n69VnyWH2m78+jtf69euxfv16dwVTqIwwmPmVeZ8R0TARLWPmXUR0DoDdOTJ2dn7vIaJbAVwJoJQw\nAOBjH4ufFHaSpyExKXzlVL2H4UsYWQeOfI2IVCht25QmDF8ZsVVoeQcjpbz6kLEeGQEWLZqqi4SD\n4du/Wc6BpDe+dm2cHEki9SGvLEflxAljo1zPMxU5Xuk2DQ0NYWho6NnXN998s7uyaC4ldRuAGzt/\n3wjgi+kLiGgOEQ12/h4A8CqYmyA6IWsgfL1OqQgja9HZELgtG/A+cvKe2eBrRPIIOaR/02PtS15S\n/ZslR8rBCCHBLMMoRV6x6ccmnB0ge4ykytMlxlrSTvXKpvd/B/BKInoCwMs6r0FEy4noy51rlgG4\nh4geBPB9AP/GzF93/QIJL0LKoOXp4pO/7e/PztVLbO7aVF3sgglJSVWV8msiJZWnS9MpqbQMKfLy\nNWpZEUbd6VQgew1IRqSxUZOUnapi07uylFQRmHk/gFdkvL8DwGs6f28GcFnod7QpwsiaFL6TfMYM\nE6am77U/MgKsWROni62ccSUvYLJ/k89s8DnkBphrt2yZ+l4IeUmMdZUVcZIbzXXrInEOA6g+wpAg\n9SZ0aZOdcsEZddI7NIzO8uolvAjfcDEvlK7biOTJkdj09t1nktIlq299bqZYpEsIYVSVHmtq0ztr\nrCX3AiVSUlL9UneEMTCQ/Ux6JQwP5HnSPh1ob6SXLk0MGdCsxRJipLNy9RKE4TuxJNJAWV5nU7oU\nbT663o/K6pKWI7XRLEVeTaSkstImEusIkIngmtpPkXBs+/omU9ZpOb2yh1E5JEI9QCb3Om+eHGFU\nsVkXGu3ELjqJvs3TRSIlFapLVuol1qDZh+G4ntC2ciT2UyRSdRIpqbPOii/DBmTWUd7BSKmUlO+8\nk3K+ytCzhCER6gEyntFZZwGHD099T8qTDvHSsiaWBHlJ5F5DyKuqlFTIfJk3L3usYw2aNSI+D8OR\nIMH+fqNLMi3L7L/PJOEczJ1rvjd9Ylyyf11hD0ZKVKFJzDsp56sMZxRhhEYYEp5RrBEBso2aRClr\nyMSSqu6I7Vsgv198I4x0HjhEF4mxlvCAgXznwKdNyWILi2PH/G6mCMiMdV/f6W2yh9yacA6y1lIT\nZbVWF4m1VIaeJQyJlAlQbYTRxB5G3iSXSEk1UWoJ5JdJ+uaBZ80yxjBGl6yxDs2PJ7360DRF7GE5\n4PSxbiq1C5zevyMjYUUSVdgGWyQhcYZII4ya0SbmliQMicqvWK+oSJfYCKOpTW/gdM8zNCV16NDU\n93zH2j6z4cSJOF2ycv4S/Ru6jqogjKbWEXB6m44f93viZJ4uUpmQkL4pwxlFGCGTIs3cNjT3qZxp\ne4QhtW/Q1B6GlJeW5Uk3kZIC2mWkJYhUIrULyBGG1EZzsn+b2n+zuiT7194dOHlmSwI9SxhVVR80\nldfO0sX3DptWhoTXKWGkLenGetJSXlrakEimpGJJMNSISOiS7hdJ8vLVJV1UILGOrC6xEZzUOpJI\n+YUUSbigZwlDMtRLexFNEUZaju/toYHJhZvMj0sYESsnNuUnRRiHD5v+8kGVKSkJ77UtKSmJ9I2V\n00SEIbHPBJzepqayBlaX2HXkgp4ljKwKkUOH/I2IxEDMmjX5pDOLkMmVNkYhC3f6dOPZJw1A6KKL\nTdUBMoSc5aWFEkYVEYaEUWs6JZXsl8OHzVz0QZtSUlURhmRRjYRjK71/AfQwYWRNipCJLpGnJMoO\nX2MneqgXIeW9Jsnr0CHTt74hcJqQQyu2kjLsMyl8DrkB1aSkQso+gdOdg9AopQrysmPtq4tEtCNB\nGFlnZUIJWUKXQ4emRvyh/ZteRxpheGD27NO9+pCBkPCAgWo8o9BJIRFKpw1aSN8Cp3ueEoflbHTh\nS15ZKSlfI5Lul9FR/7JPINtIJ2/06AKpCCNrrH2jt6oijBAHzhr6pJE+fDhMTqxtmDnTVFUl512o\nndKUVASIZIxaFnP7GhFAxtinZRw8CJx9tr8uWYQhYURCCENi0UnpkpV6aYrU0yR48KB/m7LShiGV\nM/Pmme+3iHG80ntnvmtJon9nzjSp2eSZGwlCDi1jTY+1RASnhBGAtCEJyWunI4yQtBYgF0pLGMas\nvZDYsF4qwgglr6RBCxln4PSUVAghS4yzlSPh7GR5476RV9Y68tVlxgzzY430+LjxqiUIObR/04Ts\nSxhZ6yjUOYgda4l15IIzijAkIowDB/wnFlBNSipkkgPGCCYNbKguUhFGbOpl/ny5CCOZGgjp3zlz\nTBrUPo0wdPMxi5BDdDl1arIgISbakRpru5ZsWsunwg+QI4xk/46NmUN3vtHO/Pnx6wiQWUtpxzbU\nNpThjCGMkyfN4unv95ORNmihaaA2EUZ6ooekXqSMSFqXEEKWiCSB01NSIf2bLnCQMiIhKSkio/+B\nA+b1gQNhc1diDwOYupZC524VEYaNmHwjr7TjFbKOgKn9a8nLl9jTuoSOdRnOGMIIreJZsGBywQEy\nEYb1+nwO3FkZsUYEmGpEAPP3ggV+MqQI4+yzT9fFd6LPnm3G1T63RMKgATJGLVQXqZRfsn9Dxtnq\nIjHWybUUatCqSPlJOV7798f3b2jBxtlnm++3CHVsy3DGEEao1yk1EMmJbidorEcTkqZIyzl+3BCY\nL3mlq9AkDBpznBds2yRh0IBwQ5Kcd/v3AwsX+svISlOEjrV0hBG6j7dgweRakoowQtuUNPahjleW\nsxNLGJKOl6akPJEciJh8/4EDk9UdEhHGvn1hE6u/f/JZBECcZ5T2On3Jy6Y7JI30sWOTzxnwRXKs\n9+0LM9JJg8YsY9RCvc6kLkCcUbNjtH+/DGEcOBBv1KQII3SsFy40/wuEk3FWhBFLXjGp3cOHJ2/P\nrympACQX3b59wKJF/jJmzjQGzG7WSUz00ElONHWiN2lEgKm6hMpJRnAxkzy56CQI4/hx098h5CXh\nHCxcCOzda/625CWRkgrp3zR5ha6ldIQRGqknySuUkNPrKDSSTBvpUF1sv4RGpNOmmX2PpIOshOGJ\nRYsmF93evWGTHJhq1EIHYt68yYUbOimAqZMrJhyPzWsDp/fv4sX+MiS8TmCqAZAgjBhd0impkP5d\ntGiyPSMjJvLyTRsCMnsYyXGemAiXk+zfmEj96NFJIx0z1klnJ0QXa6ST0WTIepS0U7HOQRl6mjAW\nL546EKFGWmKiJ3UJneRWFzvR9+yJN9KhBg2YatQkdImZ5IsXGx0AGcKIGaOkgxHrATOH9y0gUyV1\n1llmr+r4cUOkg4N+z3ywkFhH06YZQt6/3+y9HT0aFnklHa89e4AlS/xlADJrKUkYe/bIEIaW1QZg\n0aJJIxLD3DbPzgzs3h02uZYsMf8LhKcpgKkTPVSXxYsndYkx0sm0SWj/Jo1rTL8k2xSbMmEGhoeB\npUvDdEmOdagR6e83UcXoaBxhJCOv0LEmmnQOYh0va9B2747vX0s6vmc5gKn9EtO/1lEZGzNjFVJY\nIxVhJAk5ZpyK0POEITEQ1jAePmxuKRCSGkgbkZiU1L595qTswYNhcpYtMwYRiJtYEv27aJHpj/Fx\nYNcu4JxzwnSRiDBmzzbjOzISTsaAMYS2f2MiODvvYgxacqx37Qo30tb5CiVjwPSDnS8xuti1FOt4\nWcKIGetly0xbbForhLySkXpoaheYdJqOHjXrKYS8ytDThCGVBjrnHGDnzjiv0w4mc5yRtkZk3z4z\nQadP95cxb545BzI6atq1fHmYLpYwTp40RjYkNTBzpvm/vXvNwlu2LEwXCcKwcoaH5SKMGMNojVos\nYezaZf7euTOOkPfujXO87DoC4vt3zx4ZJxCIS0ktXWr6V2IdATL9Ozxsxl364UlAjxNGMg+8e3f4\nQCxfDuzYMTkQIejvN8bx8GFg2zZg5cowOStWANu3x3lFRJNe8I4dcRPdenoLF4Z5V4D5/p07ZQhj\nZMSkB0K9q5Ur4/vX9u2pU0ZOqJFetsz0SwxhnHOODGHYCCNmHa1YYeYbELeWrPO1fbuRGYLly83/\nW9sQG8HFrCObShofjyOM5DoKJeMy9DRhzJxpDMfu3cAzzwCrV4fJsZMrxisCJj3PGMJYuRLYujXO\niAAyhLF6tdFl+/ZwQwSY/92xI44wliwx7dm61egV6l2tWGHGRyLCGB42xiBkgxiY7N/h4biUyc6d\nk0QaEgUCZoy2b49bR0uXGoN46lR8/9qxXrUqTIbth0OH4gjDRhgx62jGDPP9O3eauRcqR2IdlaER\nwiCiNxDRI0Q0TkRXFFx3LRFtJKJNRPTukO86/3zgySfNQMQQhh2IGMJYtQr48Y/jJvqqVeb/n346\nXAYwmarYvj18gp57rmnP5s3AmjXhuthQOmair1kDbNli9AkdZ2AywoghQdueGMcAMO145pm4/j37\nbJN63LLF6BVKpOefP9m/550XJmP6dBOJPv64ceZ87+tmYddATP8SGTlPPmki5JgoMJYwANOnmzeb\n/g0da4l1VIamIowNAF4H4Jt5FxDRNAAfBnAtgEsAvIWILvb9ojVrgG9/23h6IYewgMk00BNPABdc\nECYDAC66CHjoIePthYaddrFs2mTkhWLtWtOeHTvCF8uqVeb/N20yBiUUK1aYNm3eHG7slywxpZ8P\nPxxHGDbCePxx4DnPCZMxd67xYL/3vfCUCTBJGE89ZcYrBNYwfutbcbpYwnj66XDCAIwOd90FXHhh\nuIyLLjJzN8bxAsz/3n23cXxC9gIBMy6bNsU5XoDp0+99z8ybkKIaYOo6OvfccF2K0AhhMPNGZn6i\n5LIrATzJzE8z8xiAzwN4re93rVljJmjMJD//fGNEfvAD4HnPC5dz0UXAV78alzJZtMjUw997bxxh\nPP/5wG23GQ80NN8/c6Yx1N/6VlyE8cIXAl//ukmZhBo1IkPm3/hGnBG54ALgwQeNlxYzZy6+GLj1\n1jgHY/VqY6RjCAMw/XvLLeZ3KM4/f9IDjjFGz30u8PnPx62jCy80hPHUU3G6nHcecPvtceT13Oca\nwrjvPrOmYnS56664dbR2rZm33/pW3FgXoc17GCsAbE283tZ5zwvPf74xIpdfHq7IrFnG2/ze9+IJ\n4847gZ/5mXAZfX3Ai19s2hRDGM97nom8fvInw2UApl+/8Y24xXLFFcB3vwtcemlcZcdzngN87WvA\n1VeHy7j6auCb3zTGOtTrBIwhuftu4CUvCZfxohcB3/mO+TvmENall5q5e0Vu8rcca9aYaPLo0Tgj\nfdVVZt5dckm4DLsX8tBDwGWXhct5yUvi19GcOcZpuu++uP616+gFLwiXMX26mXf33tuFhEFEdxDR\nhoyfX3AUweWXlON1rzMpgl/+5Tg5L3uZMa4xYf0rX2l+Dw3F6XLVVebE66WXhsu47DKz8F/xijhd\n3vpWIyeGBNeuNQT2trfF6fKOd5jfL31puIxFi8zCe/vb43Sx/x9DGIODhije9KY4Xa6/3vTxz/5s\nuIz+fuDNbwZe//o4Ir3mGvP7F38xXAaR6d81a0yUGwq7Hu28CcWv/7ohi5gn3L3mNeb/Y+fdG94A\nvPzlcXaqCMQsYpfDvpzoLgB/wMz3Z3x2FYB1zHxt5/V7AUww859nXMs33XTTs6+HhoYwlLDKobc2\nrwInTphJHuNJT0yYErzQ6htJMJtnErSlf0dGwp65nsTERHiJcBLM8bXwo6PGWFdRU++LkydNm3yf\nCV4Vxsbi18CJE+1pT1XP4U5i/fr1WL9+/bOvb775ZjCz8+xqA2H8ITPfl/HZdACPA3g5gB0A7gXw\nFmZ+LONabrIdCoVC0Y0gIi/CaKqs9nVEtBXAVQC+TES3d95fTkRfBgBmPgXgnQC+BuBRAF/IIguF\nQqFQ1INGIwwpaIShUCgU/uiKCEOhUCgU3QclDIVCoVA4QQlDoVAoFE5QwlAoFAqFE5QwFAqFQuEE\nJQyFQqFQOEEJQ6FQKBROUMJQKBQKhROUMBQKhULhBCUMhUKhUDhBCUOhUCgUTlDCUCgUCoUTlDAU\nCoVC4QQlDIVCoVA4QQlDoVAoFE5QwlAoFAqFE5QwFAqFQuEEJQyFQqFQOEEJQ6FQKBROUMJQKBQK\nhROUMBQKhULhBCUMhUKhUDhBCUOhUCgUTlDCUCgUCoUTlDAUCoVC4QQlDIVCoVA4QQlDoVAoFE5Q\nwlAoFAqFExohDCJ6AxE9QkTjRHRFwXVPE9HDRPQAEd1bp44KhUKhmIqmIowNAF4H4Jsl1zGAIWa+\nnJmvrF6tdmL9+vVNq1AZerltgLav29Hr7fNFI4TBzBuZ+QnHy6lSZboAvTxpe7ltgLav29Hr7fNF\n2/cwGMCdRPRDIvr1ppVRKBSKMxnTqxJMRHcAWJbx0fuY+UuOYn6amXcS0WIAdxDRRma+R05LhUKh\nULiCmLm5Lye6C8AfMPP9DtfeBOAoM/+PjM+aa4RCoVB0MZjZOe1fWYThgUxliWgOgGnMfISIBgC8\nCsDNWdf6NFihUCgUYWiqrPZ1RLQVwFUAvkxEt3feX05EX+5ctgzAPUT0IIDvA/g3Zv56E/oqFAqF\nouGUlEKhUCi6B22vkioEEV1LRBuJaBMRvbtpfWJBRJ8iomEi2pB4bwER3UFETxDR14lofpM6xoCI\nVt8yi84AAAVUSURBVBHRXZ1Dmz8iond13u/6NhLRbCL6PhE9SESPEtGfdd7v+rYlQUTTOgdpv9R5\n3TPtyzoo3GPtm09E/0xEj3Xm6E/5tq9rCYOIpgH4MIBrAVwC4C1EdHGzWkXj0zDtSeI9AO5g5osA\nfKPzulsxBuD3mPl5MOnI3+mMWde3kZmPA7iGmS8D8EIA1xDRS9ADbUvhdwE8ClPyDvRW+7IOCvdS\n+z4I4CvMfDHMHN0I3/Yxc1f+AHgxgK8mXr8HwHua1kugXecB2JB4vRHA0s7fywBsbFpHwbZ+EcAr\neq2NAOYA+AGA5/VS2wCsBHAngGsAfKnzXi+1bwuAhan3eqJ9AOYB2Jzxvlf7ujbCALACwNbE622d\n93oNS5l5uPP3MIClTSojBSI6D8DlMAUNPdFGIurrFGkMA7iLmR9Bj7Stg78G8J8BTCTe66X2ZR0U\n7pX2rQGwh4g+TUT3E9Hfd6pPvdrXzYRxxu3Ws3EDur7dRDQXwL8A+F1mPpL8rJvbyMwTbFJSKwH8\nLBFdk/q8a9tGRD8PYDczP4CcUvhubl8HP83MlwO4DiZd+jPJD7u8fdMBXAHgo8x8BYARpNJPLu3r\nZsLYDmBV4vUqmCij1zBMRMsAgIjOAbC7YX2iQEQzYMjiH5j5i523e6qNzHwIwJcBvAi907arAVxP\nRFsAfA7Ay4joH9A77QMz7+z83gPgVgBXonfatw3ANmb+Qef1P8MQyC6f9nUzYfwQwIVEdB4RzQTw\nJgC3NaxTFbgNwI2dv2+Eyft3JYiIAHwSwKPM/DeJj7q+jUS0yFaYEFE/gFcCeAA90DYAYOb3MfMq\nZl4D4M0A/p2Z34YeaR8RzSGiwc7f9qDwBvRI+5h5F4CtRHRR561XAHgEwJfg0b6uPodBRNcB+BsA\n0wB8kpn/rGGVokBEnwPwUgCLYPKJfwTg/wH4JwCrATwN4I3MfLApHWPQqRr6JoCHMRn6vhfAvejy\nNhLRCwDcAuOE9cFEUH9JRAvQ5W1Lg4heCnNLn+t7pX1EtAYmqgBM+uazzPxnvdI+ACCiSwF8AsBM\nAE8BeDuM7XRuX1cThkKhUCjqQzenpBQKhUJRI5QwFAqFQuEEJQyFQqFQOEEJQ6FQKBROUMJQKBQK\nhROUMBQKhULhBCUMhUKhUDhBCUOhSICIFnaeh/AAEe0kom2dv48Q0Ycr+L7PENFmIvqNgmte0nl+\nwYa8axSKOqAH9xSKHBDRTQCOMPP/rPA7Pg1zq/B/LbnuXJjHFL+gKl0UijJohKFQFIMAgIiGEk+Z\nW0dEtxDRNztPabuBiP6q87S224loeue6FxHR+s7tsr9qb/KW9x2d/3kDEW3oPLnv7qxrFIqmoISh\nUIRhDcyDhK4H8H9gnlr2QgDHALymc1fevwXwS8z8EzBPU/xvDnL/C4BXdW6T/guVaK5QBGJ60woo\nFF0IBnA7M48T0Y8A9DHz1zqfbYB5auJFME/cu9PcpBfTAOxwkP1tALcQ0T8BKExTKRR1QwlDoQjD\nScA8NImIxhLvT8CsKwLwCDNf7SOUmX+LiK4E8BoA9xHRi5h5v5TSCkUMNCWlUPjDZT/hcQCLiegq\nwDw4ioguKRVMtJaZ72XmmwDsgXl6n0LRCmiEoVAUgxO/s/4GTn+sJTPzGBG9HsCHiGgezFr7awCP\nFnwHAPwFEV0IQ0p3MvPDsQ1QKKSgZbUKRYPolNX+GzP/S8l158GU32pZraIxaEpKoWgWhwD815KD\nez8D86jQPbVppVBkQCMMhUKhUDhBIwyFQqFQOEEJQ6FQKBROUMJQKBQKhROUMBQKhULhBCUMhUKh\nUDjh/wM+WjlCNZ0K4gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f16e8ea0e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "solver = ode('cvode', rhseqn, old_api=False)\n",
    "time = 0.\n",
    "solver.init_step(time, initx)\n",
    "plott = []\n",
    "plotx = []\n",
    "while True:\n",
    "    time += 0.1\n",
    "    # fix roundoff error at end\n",
    "    if time > 60: time = 60\n",
    "    solution = solver.step(time)\n",
    "    if solution.errors.t:\n",
    "        print ('Error: ', solution.message, 'Error at time', solution.errors.t)\n",
    "        break\n",
    "    #we store output for plotting\n",
    "    plott.append(solution.values.t)\n",
    "    plotx.append(solution.values.y[0])\n",
    "    if time >= 60:\n",
    "        break\n",
    "plt.plot(plott,plotx)\n",
    "plt.xlabel('Time [s]')\n",
    "plt.ylabel('Position [m]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solver interpolates solutions to return the solution at the required output times:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "plott length: 600 , last computation times: [58.60000000000056, 58.700000000000564, 58.800000000000566, 58.90000000000057, 59.00000000000057, 59.10000000000057, 59.20000000000057, 59.30000000000057, 59.400000000000574, 59.500000000000576, 59.60000000000058, 59.70000000000058, 59.80000000000058, 59.90000000000058, 60.0]\n"
     ]
    }
   ],
   "source": [
    "print ('plott length:', len(plott), ', last computation times:', plott[-15:]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple Oscillator Example: Internal Solver Stepwise running\n",
    "When using the *solve* method, you solve over a period of time you decided before. With the *step* method you solve by default towards a desired output time after which you can continue solving the problem. \n",
    "\n",
    "For full control, you can also compute problems using the solver internal steps. This is not advised, as the number of return steps can be very large, **slowing down** the computation enormously. If you want this nevertheless, you can achieve it with the *one_step_compute* option. Like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEPCAYAAABRHfM8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztfXmUnVdx569aUmu3ZO27Zcu2bINXEseGSWiDcWxMzJIA\n5kwCCQkhZDhkmISDYUgsz5wJk5AZGAIkzoRtMhyWCTHY2AbbYLGEE/Ai20K2vEq29n3tbqm3mj/u\nu+6vn77lLvUt77l+5/Tp7tev69W9t279quouHzEzFAqFQqEoQk/dCigUCoWiM6CEoVAoFAonKGEo\nFAqFwglKGAqFQqFwghKGQqFQKJyghKFQKBQKJ9RKGET0BSLaQ0QbM/7eR0RHiGhD6+tjVeuoUCgU\nCoPJNX/+FwH8LYD/k/OeHzLzDRXpo1AoFIoM1JphMPOPARwqeBtVoYtCoVAo8tH0NQwG8EoiepSI\n7iKiC+pWSKFQKF6qqLskVYSHAaxk5gEiug7AtwCcW7NOCoVC8ZJEowmDmY8lfr6biD5HRPOY+WDy\nfUSkF2IpFApFAJjZuezf6JIUES0mImr9fDkAaicLC2bu2q+bb765dh20bdo+bV/3ffmi1gyDiL4K\n4NUAFhDRNgA3A5gCAMx8K4DfAvA+IhoBMADgxrp0VSgUipc6aiUMZn5Hwd8/C+CzFamjUCgUihw0\nuiSlMOjr66tbhdLQzW0DtH2djm5vny8opI7VNBARd0M7FAqFokoQEbhbFr0VCoVC0RwoYSgUCoXC\nCUoYCoVCoXCCEoZCoVAonKCEoVAoFAonKGEoFAqFwglKGAqFQqFwghKGQqFQKJyghKFQKBQKJyhh\nKBQKhcIJShgKhUKhcIIShkKhUCicoIShUCgUCicoYSgUCoXCCUoYCoVCoXCCEoZCoVAonKCEoVAo\nFAonKGEoFAqFwglKGAqFQqFwghKGQqFQKJyghKFQKBQKJyhhKBQKhcIJShgKhUKhcIIShkKhUCic\noIShUCgUCicoYSgUCoXCCbUSBhF9gYj2ENHGnPd8moieJqJHiejSKvVTKBQKxTjqzjC+CODarD8S\n0esBnM3M5wD4QwB/V5ViCoVCoZiIWgmDmX8M4FDOW24A8OXWe38GYC4RLa5CN4VCoVBMRN0ZRhGW\nA9iW+H07gBU16aJQKBQvaTSdMACA2n5n6Q949lngne8ETp6Mk/PJTwLf+EacjB07gLe/HTh2LE7O\n3/898KUvxcnYvx9461vN9xh87WvApz4VJ2NgALjxRmDbtuL35uHb3wY+8Yk4GUNDwO//PvDMM3Fy\nfvAD4CMfiZMxOgr80R8BGzNXAd3wk58AH/pQnAwA+OAHgX/91zgZP/858J73ABw502+6Cbj33jgZ\njz8OvPvdpp9j8PGPA7fdFidjyxbgt38bOHEiTs6nPw18/euB/8zMtX4BWA1gY8bf/h7AjYnfNwNY\nnPI+vvnmm1/8uv/++9kHH/0oM8B8111e/zYBJ08aGcuXh8tgZv7EJ4ycr3wlXMboqJExZUqcLrfe\nauR89rNxcqZNM3JGRsJl3HmnkfGXfxmny8UXGzmDg+EyvvtdI+Omm+J0ueIKI2f//nAZP/uZkfG+\n98Xpcv31Rs7WreEyXnjByHjb2+J0+YM/MHIeeyxcxq5dRsa118bp8oEPGDk/+lG4jBMnjIyXvzxO\nl7/4CyPnjjvCZdxzz/0M3MwzZxpfaSjA3V83PcO4HcA7AYCIrgBwmJn3pL3x5pvXYd0689XX1+f1\nIQ88AFx4IfDgg+GKbtoEnHMOcOBAXATw2GPAuefG6bJlC7BiBTBrFrB3b7wuGzaEy9i/H+jtNfo8\n/3y4nPvvBy64AHj44XAZJ08CTz5p2rRpU7icBx8EXv5y4KGHwmUwA7/4hbG7mP59+GHg/PPjdAHM\nWMfKWb8eeNnLgEceidPloYeAM86Iy5oeeMCMUYy9AKZfzjvPyAvF448Dy5cDW7cCY2PhcjZsMHMg\nxjcsXtyHNWvWYWxsHT70oXXe/1/3ttqvAvgpgLVEtI2I3k1E7yWi9wIAM98F4DkiegbArQD+OEvW\nzp3hejzzDPCmNwFPPBEu44UXjGGddRbw1FPhcjZtAt74RuP0Q7FliyGvNWvi5GzeDFx/PfDcc+Ey\nnn4aWLvWOOlnnw2X8+yz8bo8/zywbJlxJDFynngCeMtbDPmE4sABYPJk4Mor4+xl61bgDW8wYxWK\n4WFgzx7gmmvMeIXiueeAa681/RxawhkdNQ72xhvN91Bs3Qr86q+a0u7x4+FynnnGtClmHm3YAPT1\nAXPmGD8RimefBV73OtO2UDz5JHDRRcCZZ4bNx7p3Sb2DmZcxcy8zr2TmLzDzrcx8a+I972fms5n5\nYmbOjBdC68ljY2bd4Morge3bw2QAxhBWrTKRUYxRPP008Ou/HmegW7cCq1ebr1jjev3r4xz9li1G\njzVr4uRs3Qq89rVGRmhte+tWM1HOPDOuf3fsAK64Ati9OzxifP55YyurVsXZy5YtwCWXGKff3x8m\nY/t2YMkSM047dsTpcv75wPz54XL27gXmzgUuvjiOkLdtM/17xhnhc2Bw0GTIr3pVvKNfu9bYXWiW\nzWz696qr4n3DmWeajCdkjJpeknLGrl1h/7d/PzB7tskMQmUAxqBWroxzAMePAyMjwKWXGqMIdYzW\nSa9eHW6go6PGIb7ylaZfhofD5Dz/vNHjrLPiovqtW4HLLjN9cihvI3YOtmwxkyVWl507jZzTTgP2\n7QuT8cIL4w4txhlZB7BsWbiTtuS1fHlc0CTRv9u3Gz2WLYufj6tWxQUHW7aYfjnzzLgx2r0bWLo0\n3EkDpi9mzzYlzJggMFaXriGM3bvD/m/7dlNfX7bMOIJQJ20NNGbS7dplBvP00w1xxESMseR14IBJ\noWfMAObNC98ptXcvsHixiWD3pK4+FWNgwHwtWGD6J1TOzp0yzmjnTqPHihXhY22ddIwMYDybXL48\nvCxrZaxYEZdh2P5dujR+PsbIAMbn48qVcfNx2bL4LHDXLmP/MYSxY4dpy8qVRt7ISJic3bvjdFHC\naBno7NlAT0/4dlZroAsXGmcbAksYREZOaPS6f7/5/4ULwxe99+wxhgUAixaFy7G6SMggql+XY8dM\nGeq00+IdwIoVRp9QMh4ZMba2eHFchmF1iWkPINO/O3aMk86uXeEB3LZtxrkuWBA+H/ftM+1ZsMAE\nb4ODYXLsvI7pX6vLlCnG9kKzbCWMFkIjxp07zWQDzKCGRml2IBYsCHcAVgYQTxjz58fJSOoS66QX\nLDAyQjMD2x5JXUJlWHshinMAVpeYMTp40NT7e3pkdFm61PRLyIL1yIgpqc6dK0Pqs2ebPg4J4JjN\n5y9aFNe/1kkTGfuLDQQXLYrTZcEC8/OCBfHzOjSb7BrCCM0wDh4cd0YxabCVE2NYtnwDxE26AwfG\nCSOGvCR0sYa+eHE4YRw4MD5ZYp3R/PlxMpJEGhO92jGaN89EiyFOOtkvMSU/q8uUKcZRh0SvSfKK\n6d+DB02fAOFy+vvNDrRp0+ICOEs6gOmfEDmjo6Z/Fy2K8w2WvID4eb1kSfh87BrCOHgw/P+sgc6b\nFybHrjecdlqcgR46ZNYvgLjIyDqSmEhE0knbMsW+fWElhmSGsXhxHJEuWGCc4tCQWRfxRbu9hJYG\nrJOePNnYzeHD4TIAYzcxutg2hcqxWQoQTxixbUrKiJmPSScd6uwPHzb2NnlyuH9p1yW0TUNDwNGj\nRo/Q9nQNYYROlqQDOP30sAE9fNgsEPf0xEWd7YQRMulGRkwaP3fuuC4hWz8PHzYyACMntgw0daqJ\nYEMW8pPkFRrpJXWxa0Qhcg4dircXYKKzD3UAkoRh5YQ6tSRhxDjp9vkYShhWRt2EkZzTMRmGDbyA\n8DZZXXp6lDDECCPWQO1AhETSSeMKjV6TpYHeXrPL6cgRfzmHD4/rMndumIzBQUNgs2aNywmJpJMZ\nRqguzDJy2sdIgjBCyaudMGKy7FjiSRJGaN9aXWLnYzJjil0LTAYqsYQRm2HEEnJyTs+ZE1YG7RrC\nOHIkLJKWJozeXlM7PXrUX04yqp8zJ2zSJY08Rs6hQ3K6UOv6yFDCSGYYc+aEyejvByZNAqZPH5dT\nF2GMjZk2JKPgEKdWRoYRSjxJuwsdI6CcDCOUMI4cMW0Bwgnj4MGJgdexY2FbYo8cGZcTQxh2ThOF\nnYDvGsKYNi2sA6QNFIhPGYFwh5ac/DFyksYVkxkkyUsqw4glnRg5EoRx5Agwc6apa8fokuyX0Ix0\ndNToY8e67gxDYj0lKWPGDNPGkJuoJQgjWcLs6Qkn0yNHzFoXYGTEBqTAeCDng64hjFDjaq9JhxKG\ndSJWTqwzahJhhMo4enR8wgHhjjFJyKHOqH2M6swwkjJidYm13eSiLBBOPEl7mT3bZHS+kfTYmAx5\nJUtsRMbRhjjYpP1KlDCBcOJJkldoe9oJIwQvecJIOpLQyZKcuED4gLYTRohzbTcKCccYOlmOHh2P\nimLkJCdLaL+k6RJL6rZ847te1a5LaMSYdGizZpkoemjIT0ayb4HwklSyTT09hjR823TkiGmHJS+p\nAC5kDjBPbJPEnAbi5lKSMGLLzKHoKsLwdQBDQyYKsnVtKQONiQBinfSxY2ayWkhlGLFpdIycZJtC\nHX17v8RkKjY4sHbje539sWMT+yXUASTbRGTa5Gu/7f0SOgfa2xTSv+3ONVSXNEL21WVgwOzqmzLF\n/C5FGCFy2slLqiQVgq4hjJjJklyUlTDQEKMYGTFGancUhTr69okrtYZRZ0kq2aZZs8Z3X4XKAOJI\nMNmm2bP9TyMfPSpD6u1yQhxsmqMPzeDaCdlXTvs8iiGvpC4hhNxuuyHjDJxqLyG+QYq8lDASCJ0s\nScOaPTts4fz48VMN1HdAk2c5gDjCiHVGlrysnFmzTE3adxueVBko6YxCyx1pDk2if0McSVqGEeIA\n2uWE6iIxBySIRzLbaZ8DIeWxWEefpksoeUnYixJGAiGGLjH50+SEDGha5BripCUIw5aSLHn19IS1\nSWINgzm9nBTijKTLY0CzMoyQOdAuY9assDmQ1ibf/m0PvGLIS2oOWEgSRqxviKkaJOWEoKsIIza6\nmjnTRNa+5zmOHRsvJQFhEc3AgPl8i54eI7OOSLp9sgBhDkAiwxgcNGdbbDoeKkeiX5gNiSfHWiLD\nkMomQ5x9mRmGb5vS2iNBGCFRfZkZRixh2Izf10+1224IuoYwJCZLT4/Zt+17fYVESaq/fyJhAGGO\nRCK6SjOsEAeQlkqHyEi2B5DplxDSGRgwV5xMmjT+mkSGIVmSCsmyY8tagEyG0R54SWX8IQFc2hrG\n8eP+O+IkfEP7POrpMb7Cd6zTfIwvuoYwJEpSVk4s8YQ4xoEBQ1ZJSERpodlOuy4hTrrd0G1k5IO0\nMQqJPMuo98fIic0whoZMuXLq1HhdJEpS7W2SGOuZM8Mi6TJKUpMmmcPBIW1KkqBEhhEqJ21e+6Kr\nCCMk0msvvYROulij6O8/dTAlHEDolsL2SMROXh+0R2kSjj5UTlqtXirAiNUlJMBo3+EHyPRLSHtG\nRgyB2W3Gobq096910r4PLiqjJGXl+M5rqQX4tBJxSOVAM4wWpBxAiByJtDPNSUu0KSRiTCOvEF3a\nDV3CoYXKyYpefdA+zoBcVN/f71fukAx2krrYjMXnKg0bMMWSl0QGNzpqzsUk51Ko3ZVBGKElKc0w\nhFF3iUGCMNoHM8SpSZQG0nQJnXSxhJFF6iHZjoQu7Ws7EpntpEnGUftE0mWV6kLkpJFXXQHc8eNm\n3vQkPFuIvaRF9b7zemzMjGk7eUlkGFJBqS+6hjAkFr0Bfwdgd84kByLEiaSli3VF0mm6hJakJMpA\nZTijqVNNKWV4OFwGIBeo+LapTF185Uhm6u2ELNEvIbr0958q57TT/Prl+HETeCXJK3Q+pukSUpLS\nDKMFqUVvX+OyO2fs/TdAmFFIRPXMp5ZNQiaLVIbRTjwzZpiIy2cRM22MQnaItEfBRP6Rp5STTovI\nfW0mK6qXIi8fOVm6SGxwkCCvUCedtqboE9WXrYuvb9CSVAJNiq7qiuoHB815hSR5zZhhDMWnPi6R\n7TAbfZIG2tNjFjF9Ho1a1hoG4N+/aRFw0zIMiTUiXzlSpF5Wv0gFTb5loLQSplSJ2Nd2h4bM/Eue\nZwpB1xBGXSWpNKOwk8XHSUtE9WmTP6Q+LmGgJ06YA3fJMwtAWFQfW5JKOy0eIqdJGYZEZpAlx7dN\nTVrDkArg0uZAXcGkxHyUyC6ALiKMuqKrtJ0zkyebL5+dJhIZRpqBhsiRIK+semkdEePJk6YElTyz\nIKVLyJqXRJuydklJrRHVsYbRpPUUia3lab4hlDBidZHYUgt0EWFMnWpq4z7PA5CI0prkpPN08XX2\naQYaKwOoZwdOVr+EyIktSQ0PG9JoJy+J4ECyLBub2UrtQmtShmHLuzG69PYaP+Wz2SIt+NIMIxJ2\nETPW2ftOurS6NiATAUhlGL5lIAnyytrCJ+GkfduTN0ax2aSvvZTdLz4yRkdP3fYJhLUpVhdAhrzS\nZEyfbkqkPhd5SjjpNF2IZJz9SzLDIKJriWgzET1NRB9O+XsfER0hog2tr4/lyfM1rjRDlzBQoPsy\nDF8nLZVhZKXjUuRVdVSfFemF2Euso7dj1P5sZ9/AK8+5uq7jpe3wC9ElbYzsHXE+2YGEk04jdSvH\n137Tsp06MozJxW8pB0Q0CcBnAFwNYAeAB4jodmZ+ou2tP2TmG1xkhjgjiYUtqQxDYt0gTReJDCOE\ndCTWMCROwGdNFqlyR6wjCtElr1+YTyUBH11CAq/kc+QBs4bX22si++SVIVmwtxJPbvNIs2YBu3a5\n62IP7rXD9k1aQJUGqYXmNF0k5Ejp4os6M4zLATzDzFuZeRjA1wC8MeV9DqZvIBFFhJwslcowJIwi\nK3qNdbBSJak6dJEir6xI2sfmsnSRGGvrpF13xOURhkSZzWess4KdEPIqy0lL+JdQXSRKUp2+hrEc\nwLbE79tbryXBAF5JRI8S0V1EdEGeQJ/JOzZmop9p0ya+LrGwBYSVcCScdJZjjC0n1VWSKrNfJEpb\nU6aY2rjrI2PLJHXAr4Qjle1IELJEe4rkuOqS5RukAjgJ31BXhlFbSQqGDIrwMICVzDxARNcB+BaA\nc9PeuG7dOuzYAfzDPwDMfejr68sVPDhoDKKnjTJ9CaO//9RdPICMM2pShlFnSSqvPh5TepHQhWjc\nZtLswEeXHTvidAH8gqYsGfZEfqwuTSIMn7k0OGjKaO22FTIfFy6M0yXrhHZohrF+/XqsX7/e/R/b\nUCdh7ACwMvH7Spgs40Uw87HEz3cT0eeIaB4zH2wXtm7dOmzcCLz5zUABVwDInyy+aefSpae+Xsca\nhlQKLHHSW3LRu71NyWuvXdLsvEXvvXvjdAHG+zeGMKQI2cd+8+ZA1dmkBAHmyWmSLj7zcXjYEFf7\nCe3QYLKvb2Iwfcstt7gLQb0lqQcBnENEq4moF8DbAdyefAMRLSYyPE9ElwOgNLKwkIiu6jAKKye2\nDFRmJN3ba6Id13MuUjuTynYAEoTs66TL7BcJXaRq9T5tkgzgJMgrq1+kMn6ffpFYk5HaVltbhsHM\nI0T0fgDfAzAJwOeZ+Qkiem/r77cC+C0A7yOiEQADAG7Mk+kzoFLRlYRR2HuX2neTTJ9uXh8dPfWK\njSxdlrevArV08bnZMs1Ik/vHe3uLZZS5SwoYd/aLFlWrS1lRvWT5USLDqDqql8iY8nTxddJlj5Gr\nnDI3SYSgzpIUmPluAHe3vXZr4ufPAvisqzypdLxqA826dym5f9xlO6BEfdxe1562FdI6gNNPL5bT\n33/qVkvAL6q3p2HTLkyrY3G3TGfUNPLyXceL3SXVCeQlWX2QyI591vH6+4F589w+Mw9dc9IbkImu\nrHNyPbpfphORkuMjY2jIEFeak/aNjGLLHXlbAWfO9Ns+GuvQRkdN37TvnAGalWE0KWiqQxeJiLxJ\n1Yc8P9XT41ci7vRtteKQKElZOVU6gLz6otROE5/oSkqXWMLIG6Pp06vdDZS1c8bKiSVBqVJdXYTR\nJF3KWvOaPt1cZOl6xUiZaxiScnzQdYQh4YyqjuqLdJFKXyV0ia29+jjpskldShcpe/Ep1Y2NpWeB\ndazjlTnWdh3P9YqRMp10cgt1jC4SmwEAmfnoi64iDCnH6GMUWVs7fU+55hmFRLYTmxkA/pF0mRmG\nRFRfhy4S6wbW5srKdnxkjI2ZqLusUt2kSeNXjMTIkSyFVl19kApUSt8lRUR3OMg4yMzvilclHr4G\nKuEYy047q3ZGZRtoXWMUq0vemQ+JqF7CEVldfOQsXnzq6/Z55yMjp97t1A5bqms/AGt1OXzYXZei\nNhXdSZX3GFLJ/q1yS6yULlIZRtEuqfMA/AHS73Pi1uvOu5jKRtMixjqiiDIjaamSlNQYxZY7JPsl\nlrx6e92ddJ69SJCXLb0MDhbvzsvTxXedKWsXjx2ntF13SQwPZz+GVGp9s2nl6kZlGAA+xsw/zHsD\nEf2XeDVk0KSdSb51SqmFrTIX4CVKUj5OpMgxxjp7iYVzKV18rhgpst19++J0AcZ1KSKMsvvFR04V\nurgGBzbbScuK6lj0rmQNg5m/XiTA5T1VoWoDzTpwB8gNptTibh0ZRlPWMLLkTJlixtBlC7VkmSK2\nf5vmpMsspwLuTrqKaNxVzsmT6WerrAyp9U0Ju/OB06I3Ef0yEd3WeojRxtbXY/EfL4uqJ8uJE6be\nm1W/lajVd2qGkdUmn10veWMkkakkSy8xulS9/tAkwig72KlDF4mSXxX2IlEK9YXrSe+vAPgzAL8A\nMBb/seWg6sUkqYkrYeh5u1XsIyrHxtLJLYkiAz10qFgXILtNdtdLlq7tuuSR1549brq4RMExZSCp\nsZbKMKQ2W0joUuUW6iaVpCQdfdqNt1ZO0xa9LfYx8+3Fb6sXnWwUsRlG1nXtgHlt2jRDGkVGUzTp\nXK4YGRsr3lU0MOBGGE2JXqsoSXViGUgyIy1zjKpe36yC1F11scGky5MPi+BKGLcQ0ecB3AfAHkZn\nZv6XeBXk0CQD9dn1MjCQ/rQxq4tLVF9Uo7RGWkQYeQ5AgrysLgMDxXfbFGVeTYleJUsMnbiG0Y26\nxM6BpvVL1i0FvnAljHcBWNt6f7Ik1SjCqPrgXp4Mn10v/f3p++EBGQMFqo3SJHXJ2krpU6rLugPK\nV5cmZRhNckZlrjNJ6dKkdaZp00zE71IiLiphupxzkVq/ANwJ45cAnMfsekC/HvT2mntehofT92In\nUWQULlsTXR1j7DZJCSftKifvVkufxcc8A/WZdCtXpv9NKrpydWrdWh+XWGhuUkkqSxc7zi63u5a9\nztTTM775o8iRF/XLzp3FukitXwDuV4P8FMAFMh9ZHqR2vUhEEYBfHbhsXarMMIoIQ8JJS8gAZOzF\ndZxtiXLq1GxdXPtXgrzKtrsmZTuTJ5uvkyfj5EiQOuA+1mUTqS9cM4wrATxCRFsA2C5nZr5IRg05\nWKcWG9VXbRRl1kx95cRGjEURjY+TjtWlCiJ1Hee8O6AAP7ubMydbF9dS3YkT2QuhnVqrd5kDVWy2\nkLS7suejD1wJ41qZjysfEp1YtVF0ii4SpOOjS9MmbhW6uPZv2nPkfXTJO0PkI6e/P5u8bIm4aONH\n3qloq4vLs9ddApW8dTGLJgUqUmdCKs0wmHmrzMeVjyojRqmSVFUZRlX9IrWG0SQilchIXR1ajC62\nNl60oOrSL0ePuumSRV6ud1LlPbTL6tJtJeIqM7giXXyQu4ZBRA8XCXB5T5Vo0kBIRfVVRq9F22qb\n4qRdy1qSC/B5Trq/v/j0ehVOpKfHZA5F14FXUTJxlVOFvVg5VZ2haFJmW2VJ6nwi2ljwnoyEtB64\ndGLeNchAs9Yw6sh2XkolKYnF8ylTTIQ8NJS9oO2ii7STzvusKoIdwK1/Je0la1edqxy7tlPmNmxf\nORLlsapKUuc7yBiRUEQKLk5teNiky2WnwFWvG1RBXgMDxVsTq5gsrtskXXSJ3SUFjI9TDGFUec6l\nyIlUmU1WRaQu/Ws3AsSu7QwMAAsW5OvSdRlGJ61dWLg4AEnnumhR9t8l1jCaRF7JrYl5O03yFjB9\ndMlr0+TJclH9wYNxulg5EnZXVa2+igDDymlSSapJuhT171DrTo2swFbCv/iiqx7RCsgZhc82yVg5\neU562jRjOEUPnm9a9FpUppCI6iUi6aqdUbdF9RIL+U0ijCo2JrjKqWqMfNCVhFHkGKsqU7gM6PCw\nqZv29qb/3fXB85LGFZvWN8kZdZoDaAqpS/ZLkzKvWNvtRFLXDCMHEgMhZRQ+pYG8OnxVhp53RbqF\ni1NzybxciHR0NJtIrZymOCOJQEVye25VpN6kklQnZTuuviGvb5N3UuWh8gyDiH6TiJ4moqNEdKz1\n5bBLu3pIlAZszbDoSWxVRBFAdTVpu/aQt39fSpciGUWnol3lNKnc0UkZhk8JM3a9qqp+kXDSklfS\nFLWpyNETuZV368gw/hrADcx8GjPPbn0VXL5RDyQcGlCdobsMZlXRqyt5VRG9uvZLFespNnDIu9BS\nojxW5XZuyetbivq3SaQe2y9TphhHbRekY3SR8lOx89EHroSxm5mfkPnIciE5EC5yYqMrl3RRyklL\nGFZV/eKiS1V7/JsUYLiWDZtUquumdSYpOVUGtkXrkj5wJYwHiejrRPSOVnnqN4noLbEfTkTXEtHm\nVrnrwxnv+XTr748S0aVFMjvRKKQi6dgdOC6GVVWG4UKkVTmjKgnDpVSX93CqKnUZHc1/1oiULsnH\nFhTJKXvDBtC86oOE/brC9fLBOQAGAVzT9nrwA5SIaBKAzwC4GsAOAA8Q0e3JTIaIXg/gbGY+h4h+\nBcDfAbgiT27VkXQVjrETM4zYckeTxqhqXfIOI7o4EddAJesOKMDtIT/2kFvROlPsWCfvpMorC3bi\nWFehS+W31TLz78p83ARcDuAZeziQiL4G4I0AkqWvGwB8uaXDz4hoLhEtZuY9WUKbVBuUiiKqWmiW\nyjCKdkmE9awVAAAgAElEQVRJ1LWB6kovkk4kr38nTTLRdN6141U5NCJDGnkP+XHVpeipcAMDxTfI\n2jZlPbbAPmskb1dd1Rl/VcGkRJtc4bpLaiUR3UZE+1pf3ySiFZGfvRzAtsTv21uvFb0n93ObFEW4\nloGqIK9OXPTutkjPtX/z5FSpS5GcqnXJszvXXXWxi95WTqetYVRdkvoigK8AeFvr93/feu11EZ/t\n+rjXdhNI/b9169YBAF54Adi1qw9AX6bApjmjqjKMKmumEoveVaztSJXHZs4E9u+Pl2PHKSvirirb\ncZFTNWHkjVORzUnrkifHru3kXVlTdYnYjvX69euxfv36/H/IgSthLGTmLyZ+/xIRfTD4Uw12AEje\nLbkSJoPIe8+K1munwBLGI48Av/u7+R88MADMnZv/Himj6KQMwzUFPnYsThepieta2spzJFOmmDp9\n3nPgm5RhVBUBu8hxcdISO9lcdZHasJF3aaCVk6eL5BkiiaApaTN9fX3o6+t78W+33HJL/j+3wXWX\n1AEi+h0imkREk4notwEUxFOFeBDAOUS0moh6AbwdwO1t77kdwDsBgIiuAHA4b/0CqG7iSqXAkoae\n16akYwyVAcj0rz37kPf8iKZFr1WVBopspqr1NyunqgyjiHiqLNXFnouS8C+ATHBQ9CgHX7gSxrth\nylG7AewC8FYAvxfzwcw8AuD9AL4H4HEAX2fmJ4jovUT03tZ77gLwHBE9A+BWAH9cJLdJ9duqjMLl\nYBmRW/QaS15Fj9sExh/yc/Jk9nteqrX6l+oaRqeVxyR0kSL1PDnDw2a+5fkGH/g8ovU3ZD5ygty7\nAdzd9tqtbb+/30dmk6Ir13LHkiX5coqctGsEYfsm6xnMUgaa97jNpJyBgew9/J3oAKR0kcgwmtQv\nVWyhdiXjKha9mzRGkof2gALCIKIPM/NfEdHfpvyZmfkDcqrIQLI0EGsUyZtms5y0xNY5H8IoklO0\ntiOty7x52XKytlAmZeQ5I2bz99jFUJeSSTdmGE3SpZOctIuM5GMLJk3KllNF5uWDogzj8db3hzBx\ndxLBfZdTpZg61US5IyPmITtpcDWKvP3jvo6xzKjeVReXktSyZcUyinQpcq5A8WKoS+ZVNFlOnDD2\nkHcq2sqpYg1Dqj4utRAaS4JF521cdalqDkyfXnwYsSrCSB5GnDUrTpcDB7L/XmmGwcx3tH4cYOZv\nJP9GRG+TU0MOyYGYPTv9Pa4DsXNn9t8lnXTstlrfklSMHKmIpoqIUVKXLFtylVH0rGgLibGu6oCa\nq5NuSlRvb3cdGMh30i79kveURt/5mKVLVZUQH7guen/E8bVGoIpU2jWSrmKBrGryklpPKdsZueoi\nIccl2ym6A8pFjoQurjtnuq0M5CKnSicttS4T2x4fFK1hXAfg9QCWE9GnMX6IbjaAgqvA6kNV0atL\nqidRYpBe9I6RI1GmcJHTJCdSdYBRNEZZ6z6uurjunKki2AHcnXTe+R+foKmTsuwqzqf4oCjD2Amz\nfnGi9d1+3Q7g1+XUkEWnRfVVLD666tJJGYbUGFVBpK61ZAl7sRlT1jkXyTEqIsHe3vzzPyMj5m95\nB2BddakyqpfyDVXYXWUZBjM/CuBRIvoKMzc2o2iHhAOocmdSFQeFpORURaSuusQuVkvJkXRoe3KO\nprrImTLFZBDDw+mX8fnocuhQnC7J8z9pGz9cDsAC1QVNVa3tuMiRKo9Vua32/zHzWwE8TKeOKDPz\nRXKqyKFJqZ6EUbgYuUu5o4oFVcnSSydFelOn5m+TdLUX6UAlljAks8k0wpAMvIqu9LBystrE7LYx\noUnZjtQYuaJoW+2ftL6LH9orE3kDMTpa/MSyIhmAe6onsYZRJMO13CFRNrEOaGgo3hl10xpG0TbJ\nOpx0f3/6uRof25WK6rPaVOUYFckZHHTfhl1VtlNVKdQVuV3DzHZj6T4A21onvqcCuAgZlwA2AUVG\nUfTQlyIZQLVrGNOnm8hnbCxOF+mIMUaGVP12cDC+X6pYl3F10hIlkyJdpJy0xAaHqgmjaeSVpYtd\n88l7voekLq5w3Vb7YwBTiWg5zN1PvwPgS3JqyCJv0kmVkiSMy55ELpLT0zP+MJssXaQWVF3lZBl6\nlbuk7J1UJ06Ey3DRRSKDaxqpVzkHiuZjk/pFynabQuqVZhgJEDMPAHgLgM+11jVeLqeGLCQivSqi\nCNcUGMh30pLlsU7KMIrkNCnDqLJMIaVLFXOgSQGcpC6xZ7Sq9lOucCUMENGVMA9OutP3f6tGk6II\niehKSk6eDNcU2MrJcwCxk8X1DqgiOU3LMJrkGJuybtCkkpSUk5bYQu3bL1lbqOvKMP4jzMnu25h5\nExGtAXC/nBqykIhopk0bv3cmDRLG5TOYUlFa7GQByo/qXe+AqkIXHzlVZF4SayGS2Y7r7rwmrRtU\nEWBUlXnZ58BnPSqglgyDmX/IzDcA+BwRzWLmZ5t4U62FRKpn753JWzeoykCB/MlbpRMByp90TdLF\n7qrrtBJDk6L6JukiQV72gtPR0Ww5sdmky+WQrnIqJwwiupCINgDYBOBxInqIiDpyDaMOZ5QX1Utl\nGFVmO1IRY+zCudWlzKjedVcdUL69uDycykUXqW3YEplK09Z2XPol+diCNEiVpHx8g4QcF7iWpP4B\nwH9i5lXMvArAn7ZeaySkOrBsQ5dcw4jd3SRFpK7OXjLbaQqRSthdXr+cPGlOcWdd2++jS5MWvV1k\nTJsmt7W87Ki+E4NJV7gSxgxmfnHNgpnXAxDkLVlIlG+AZjmjshfrfIm0zAxDkkirKlNI6dLba+5X\nGhmJ00XCoeXdSTU66nYHlJQudmt57BbqonkkERxIZHBV6+IKV8LYQkR/TkSriehMIvoYgOfk1JBF\nFSWpKvfmWzllOmlJIo29pqRJGYZUacBn7Syrb3x1ie2XyZPNomraxYGud0ABMsEOkN2/VZfqiuRI\nBZNVBgeucCWM3wOwCMC/APgmgIUA3i2nhiykoleJgSiKXqXKY522xVfKyMvOMOqYuFn2W7Xt5smp\nS5e0fhkaMuTmUqorqj5UGdVXsYYhnWEUXT44HcAfATgbwGMw6xiNv7W2W6OIMg/uVb3oLbmeUuaZ\nEMnFx1gS9O2XbdvidbH9234nVdWkbnWR6Je8OZD1KGVXXcbGDIEV3VXnokvVdueCogzjywBeAWAj\ngOsA/I3cR5eHsqMryYNlVS7Az5hh9E6rSVcdMTYpw+jtNX2SVnqRzDBiHUDV5dQiXXwWiMvM+JvY\nL7G76nwDuEZkGADOZ+YLAYCI/hHAA3IfXR6KJu7SpXFyfA6WSWYYWQ97dzWK5N1L7ZNdKsNw3SVl\nz7iMjZ3aj74O4MiR9L+FOJL2CLPq9RQgu3/rKo81qSQVq4vkonfsGBU5+qInKyZ1ic12XFHk8l7c\nq8HMKfs2momy65RVn4q2ctJ0GR420XHR4zaL9KnaAeTteqmr3JHlAKpc9M6T06lloLK3c0varkSG\nIbExQcJefM4QuaIow7iIiJJP0p2e+J2Z+TQ5VeRQRUTjahTWkQ8Pn+rU+/uBJUvc5BQ5elejKNMB\n+OxWsbqkTYy6ovqmRNJNIfU8OXX1S1NKUlm6NKmEKX3KGyh+RGvKs8OajyZNlqSc9nLH8ePpD9vJ\nkhEbAefJOX4cWL48TsbwsCEun2ynSZF0mRlGpzrpvOhVKsOocm2n7HMYPu2ZOnX8zE37Di+J8piP\nf3FFY2+cjUHeqVCJKMKXMPJq0nXokhWNSOjiml0A2YYu4dBGRtwPlgHVRPWxDkDCuVo53URePjKm\nTCl3g4OPjORTGtPkSGQYkgveQJcSRt4DhyQM1DfVy5IjkWH4OJE8OT665BGgT0QjlWGk6eJzsAyQ\nmXR5i4+uu+qAaspjVQYqTVpPsQcjy8ompeajBHkdP94lhEFE84joXiJ6iojuIaKUJw8DRLSViB4j\nog1E9HOfzygzMgotSbWj6pppkZzYyeKbApeZYUhlgRK6+Oyqs3IkdClzQTVkrSoNEm2qYz6WHUxK\nrWF0S0nqJgD3MvO5AL7f+j0NDKCPmS9l5st9PkDKKMo0UImovkm6hBBGWRlGk5yIVNlQgryGh03G\n47rOJOGMrIyyzv80aaxD1hSboosL6iKMG2AOBaL1/U057w3aFFbmpPNNO6VSYImSVJm6SJakYm8N\nbVKGEVKmKGsHjlSpzic4sA/5ad9CHbKrLqtffNbOyp6PVa9vSmX8LqiLMBYz857Wz3sALM54HwO4\nj4geJKL3+HxAJ5SkfKP6pqynZOniWzMtcz3l+HFg9mw/XZqUYcQ6I/uQn/ZbbyVtN5YET540RDLJ\ncS9m3hjFjrVdZ5JY36w6w5AoM7vC4bquMBDRvQDSThn85+QvzMxElPFEWryKmXcR0UIA9xLRZmb+\ncdob161b9+LPfX19mDGjr9TaYNWTrkkRTdZ2QN+IJo94XB1AVt8eO+avS2z/Jq8DT0bwknXtkB04\nyb4MyXbS7qQKHev588dfk5xHy5b565LE4KDZKOOzzlR2eayMRe/169dj/fr17gq2oTTCYObXZf2N\niPYQ0RJm3k1ESwHszZCxq/V9HxHdBuByAIWEAaQ7WOvkenvd2lCFUYTUgZPOSKrc4eMAkk8cOy1x\ndFNq0dvH2Ze9nuJDpPbG1KGhidt5pUj92LGwSDr5P8eP+9uuVP+2y5EKMHyDgzRdJDODhQvd5aS1\nyd7eEOun0vxLX18f+vr6Xvz9lltucVcW9ZWkbgfwrtbP7wLwrfY3ENEMIprd+nkmgGtgLkF0Qlon\n2okbW7/1nbhpDmBoyKTBrkaRVQcOidIkDgCmtUlqDcPHkSRP0ofKsLqUtXguRV4SDlaqVCfRplAC\nLEOXEPuXCiazyEtinalbFr3/O4DXEdFTAF7T+h1EtIyI7my9ZwmAHxPRIwB+BuA7zHyP6wdIGIVU\nuSNLl1mz/O55STNSiW21vvVbINsxVr2GkSWnjjWMLDlS5CXhYOty0mltChnn48fL0UWqPCZhL1K3\nN5Sxrba0klQemPkggKtTXt8J4PrWz88BuCT0MyQmi0SNPUuXEPa3hp6sA/tc6ZGli92p4lq/TeqS\nRIgD2LVr4mvM4ZF08pkNdWUYaXLqjKQlykBllfx8+2X2bBnCKDOYlMhI61gjdUVXnvQG0p19yMSV\niPSyMgNf9k/TR6I8FrL9rixndPKkKb+5nhOwupRVBvKdvLNmmTFJoi67mzXrVAcb4qTb2wPUk2Hk\n6RJLyJJbnyVKu3Xo4oKuJYysiMbXidhnNiQhFemFZBhlOMYQw5IgQQlHb3Upg7zsOYFYpya5eyxW\nl7oIQyLYydJFYkecVBlIYruxL3ll3Y/VTecwSoeEgfb0mIGLjdLy1jB80KQMI0uOj6GnRcBS2Y7U\nOtO0ae7nBAAZJ53WntFRs+HBJ9uZPRs4enTiayGk0y4jVE7sWEuSV2zQNH16+gWnx45N3DlYhS52\n12IsCbqgawkjyxn5TFxAJmKUjKTLWMQMddKxcrKcq68up50WP0ZppaSjR/0mPyBDGHmbG3zWmdL6\nRSKqZ/Z3RqeddirxhPTLiROGPJOQWvT2aU9PT/ocOHpUZqx9Hb0EIbugawlDYuJaORKGLsH+ZZVe\n6tIlzYmEGHmWHJ8xmjNHzl5iycs+vjZ595KULr6EPGuWGeekLoODZot3+zMc8pBF6j5tsjfNJh1j\nyCYJqbJslm/wzTDSsmPfQCWtfzXD8IAUYWRFaXVkGGVtk5TSRSLqDMkCJYhHIgK2cmLtLu3MTcgY\nZZWkfHWZPn2iU5PSJSSbbJ/XQ0Mm2nc9zwTI7VqUsBmJLAXIHmslDEdIOQCJSSe50FzGeYNuzTBi\ndQktSUmQYPs41UVewKlOuq4sME2XENJJs92jR099ImYR2vt3ZMTs8ovdEhtid2n9W8Y5jK4lDMmS\nVFKOvdBt2jR3GRK7mwCZxd0mbauVyN4AmcXdqVPNAubJkxN1qaMMBJw6TqFRfWx5zMpJ9m8oYUiN\ndSx5STnp9n6x7Yk9jCtBGCHrTC5QwvCUYw3U1yjaF6QktrKOjflv+yxzW63v5J02zZDv0FC4DEAm\nwyA61amF1JKl7E4iw5BYf7NyYp20VOYloUuWkw7JMNoJw9de0nxDKGEk++XkSWPTPqU6F3QtYUjW\npNuNIkTGkSMTXwupL7ZH9SEntHt7zS6T5LXXUml9yHZAqXKH1HpVcqxDa8mxi7tAeRlGXVF91nys\nK8Not90jR+KddIi9zJlzqm+Qsl1fAnRB1xJG2RmGD+xgJneaHDniP6DtkVHIhEveNBujS7uMoSG/\nGzYt2g29rkXvNDl1lqTKWsMI7d8ySF1qDSM2ewPkSlIShCGhS8icdkHXE0bSSddloHYL4uDg+GtH\njky8+8gFaVGnry5pciR0sQbqU6oD0p10rKO3Ja5Y8pI6h9G0DEMiaJLaUVRHhpF2RquukpQl0uQB\nQIk1jJCMyQVdSxj2SV7JrYkSC6qhh2HaIwmJqD5k8ls5ac4+Rsbhw/6kA8hFr+1jFLLYV0aGYTOv\n5PMxXNDev5I7/GKddEiA0aQ1DFsibs/46yhJTZpkxrpdjoQummF4QiI7SFsIDXHSEoTRHnWGEKCk\nLrFlLUBm3SAtugolr9iIUWLbp9VFYqG5fdvn8LDfDj8rJ9kvIcGBPQCYjKQlMowQe7FP1ktm/BJR\nvZRv0AyjJkg4ewkDtbrE1hjbo85Dh4DTT/fX5fTTzaS3OHw4nrxiMox2Z+TbpnYZhw7J6CJxgCq0\nbDh37sQxinGuNpIO2faZlGMRMtY2krY2MzpqHHbI6eqkLqFzYO7cU5207xxIW8MIcdLthKGL3jVB\nYlEqLboKMdAyMoyYyXLo0ERdfB2AZEkq1tlLkE6aHIkAI4SMgVMJI2Ss29fOQiJXQIYwrBzbvzYC\n9tnhl6ZL6ByYM2e8f+3ZG9+yoURJyuoSm2GkLXprhuGJ5EAMD5v1jNjdHTHRayxhtGcGUhmGhC4S\nTsTK8W2TlWEjaakMIyRinD594tmSmDFKknqonKSDPXgQmDfPX0b7HJBYrwptj1RwkMwwQp2rVEkq\nqcvwsP9pcauLrmFEYt688UlnjTx2F09MRGONwi7E+9aS580zkz5Wl2SGwRxOGAcPjjvpmHWD2Ihx\nyhTzZSNpqQwjJGIkmpgdxIyRRHCQbFOMjGSwI5FNShJpiC7JDCMm85LYyprUxQYpsX5KM4wAWKcG\nhEdXSRlAXEkqmY6HGJadLMlIOjbDGBw09WXfdHzaNOOkbVkqtPSSNHRmOWckVR4L6d/588dtpm7C\nSDrYUBntgYrEGIXOx2TfWl1iM4zQaLy9lHTw4MTHJ4fIiVkjjV2TcUFXE0bS0GMnv0S5I5kChwzm\nlCmm5BGb1iczjJgDPsn+jXEitl/6+03dPeQ6g6QciWicGThwIMwBSNidFGHMn2/aAYQ7aSnCSNpd\nDHnZ9lg5dWUY7f0S2r9Jwgi1OWsvST+lhOGJZHQVOpi9vSaajnXS7c5VwklLRJ1ShBFaklqwANi/\n3/wc2p52XWIcmpXR328WjH3Lhu26SJD62Fh4/yYJQ4J0gPCofuHC+LFu10Vil1Qo6Zx+upFhH+hU\nJ2H09ppg0gY8Bw6Y/pZGVxOGxMS1cqyRhjqjhQuBffvMzwcOGEcZq4tE9CpFXhL9EjNGCxcCe/fG\nyVm0aOIYhUx+QD7DOHbMOIMpU/zlSBNG6JoXYGze9m+oc505c3wDCyATHOzbF+ZcJ00ymUlsUCpB\nGMDE4Gv//nAfk4euJ4zYwQRkatJJh7Z3r3FOIZDOMPbvD49E2nUJcSJJwgid/MBEZy9BXqH1aEB+\njGKIVIIwZswYPzdx/LjJunyetmchkU0Sjc/H0VFDpiF2l7SXGOcqUfJLyoghjKQcJYwAJBesJTIM\nW9cOMYqkM6qbMJLljr17ZQhj376wNklmGLFy5s83/zs6GjdxJcZo9mwTSQ8OyhFGqENLOul9+8Id\nkdRY2/49csRsk580yV/GokXjAVxohgGM9+/wsLn5IGQtJKmLRIbBbL6HyslDVxNGsnwTU2Kwk+XY\nMXPQKGQXQzKikSAM5rhsx+oSO1msY9yzJ6xNc+YYp3jypJwzCu3fyZONPgcPypWkdu8GFi/2l0EE\nLFli/j+mX9qj15g5cOCA6duQ9gATM4yDB+NJMDRIASY6aYkMw85F3+2wgOlPCcKwutj1t+nTw+Tk\noasJY+lSYOdO8/POncCyZWFy5s83RrVrl5nEIUiyfyxhHDhgyi7Tp/tvh7Uyjh8fd9Kxupw8aeTF\nlBj27zfOcenSMF0sYTAbOaHjZB1JzMRdtMjoABibCW2TJYwdO4Dly8NkJAkjRo4d6z174gjDkvqe\nPeFjZHWJsZekk5Yg5NDAADD2smeP+TmWMPbvL68cBXQ5YSxePJ4uxkyW5cvN/8cY6NSpphZ8+HAc\nYVgnEuOIenpM3+zeHVeSsjJsluJ7zYOFdfYxhGzXiI4fNyQU+ixjq0tMgLFihbGXkRETSceOdYwu\nSWcUMwds8BVDGMldUjt3htuvJXWJwACIW8ez/RJL6kePGj8VM6/tHCirHAV0OWFMnjwe7e3YET7p\nVq4EXnghzqEB484oJq1fuRLYvj3OsADzv7t3x0VGK1cC27bFORFgImHEOJF9++KciJWzd68Z71Wr\nwmSsWGHGaM8eM3FDFoiBiYQR6oxWrwaef944pLGx8NO/Z5xh5ISWHoHxDGNsLG6sV6wwdhcz1jNn\njj/3eseOcF1WrTK6bN9u9ApBT894dvDCC2ZeheCMM8z/x5BXEWohDCJ6KxFtIqJRIros533XEtFm\nInqaiD4c8lnLl4872BjCsAYa46SXLRt3JKERjZ0sO3fGOcYlS0yfbN1qnEqMLtu2xRnosmXGyGMI\neckS0yexpG6zphjCWLrUjPH27XH2IpFhzJ1rHOOmTWaMQmrswDhhbN8ePtbTpxvC2rzZ/Ox7X5KF\nnY8xY23XiJ57zqw/xAaTsU562bJx4okhjK1b42y3CHVlGBsBvBnAj7LeQESTAHwGwLUALgDwDiI6\n3/eDli0DHn3ULFSHHMICxg10+/ZwwwKAs882E3f//jgD3b4dePZZ4KyzwnVZvRp4+mljXDGEsWMH\n8NRTwDnnhOty9tnAM8/EtWnZMhNFb9gArFkTp4vtl9BJ19trIsYNG+IDjG3b4pwRkRnfn/40zqFZ\nwnjmmbixPuss4P77w6NxYHw+Pv103FivXm10OeOM8HKq1SWWMM49F/jJT0wpNZRI7Rh1HWEw82Zm\nfqrgbZcDeIaZtzLzMICvAXij72edeSZwzz3meyhWrDBR3uOPm4ENxdlnA/feawY2tEyxaJHZvvfA\nA3ET92UvA+67zywghu6mmDHDLHT/4AembaE45xzg3/7NLJ7HrMucfTbwne8A550Xrsu55wJPPGEc\nQIxTW7UK+O534+xl7VrgyScNkZ5xRricpGMMRZIwYsZ6zRrg9tuN/cXIeOop0zcxY71mTbxvWLXK\n9Muzz8bJWbsW+N734hz98uWmnPrEE+FBYBGavIaxHMC2xO/bW6954aKLgNtuM99DMW2aMYa77wZe\n/vJwOeedB9xxR9xk6ekBfvmXgbvuipNzwQVmspzvnbNNxBVXGENfuzZcxtq1wPe/b/ontGRi5dx3\nn4wuZ5wRnpECwMUXA9/+dpy9XHCByQyY48psa9cCd94JXHJJuIw1a0wpKabGDhh7u+ce4MILw2Ws\nXm3WHjZtiguaLrnE9EsMqc+fbwKn+++Pa9PataZfYvzU5MmmP+68M05OHkojDCK6l4g2pnz9hqMI\nLn5LMa65xnx/05vi5LzmNeZ7jIG++tXm+xveEKfLVVeZ75deGi7jFa8w31//+jhdbL9ceWW4jIsv\nNt9f+9o4Xa6/3nz/tV8Ll3HWWcYBxMgAxsfY2l8IFi40zvnqq+OI9IYbzPfrrguXMX06cNllwCtf\nGZ4dA+P2dvXV4TKITKBy7rnh5RsAuPbaiTqF4uqrzdpXTPnRzqNY32DbEpN55YGYRfxy2IcT3Q/g\nT5n54ZS/XQFgHTNf2/r9IwDGmPmvUt7LN99884u/9/X1oa+v78XfBwbiDAsYv78m5NBeEnbLW4wD\nGB012zVjLxc7dMgsQoaclLUYGzNyYrfxHTxoDs3F6GJP4sfuQT90yOyiCbk1N4mYix0tQp8GV5Yu\no6Pxcynmmh6LwUGjS+j2aQuJMwtDQ0aX2INy9mrzGN8wNma2lmfthlu/fj3Wr1//4u+33HILmNn5\nE5tAGH/GzA+l/G0ygCcBvBbATgA/B/AOZn4i5b1cZzsUCoWiE0FEXoRR17baNxPRNgBXALiTiO5u\nvb6MiO4EAGYeAfB+AN8D8DiAr6eRhUKhUCiqQa0ZhhQ0w1AoFAp/dESGoVAoFIrOgxKGQqFQKJyg\nhKFQKBQKJyhhKBQKhcIJShgKhUKhcIIShkKhUCicoIShUCgUCicoYSgUCoXCCUoYCoVCoXCCEoZC\noVAonKCEoVAoFAonKGEoFAqFwglKGAqFQqFwghKGQqFQKJyghKFQKBQKJyhhKBQKhcIJShgKhUKh\ncIIShkKhUCicoIShUCgUCicoYSgUCoXCCUoYCoVCoXCCEoZCoVAonKCEoVAoFAonKGEoFAqFwglK\nGAqFQqFwghKGQqFQKJyghKFQKBQKJyhhKBQKhcIJtRAGEb2ViDYR0SgRXZbzvq1E9BgRbSCin1ep\no0KhUCgmoq4MYyOANwP4UcH7GEAfM1/KzJeXr1YzsX79+rpVKA3d3DZA29fp6Pb2+aIWwmDmzcz8\nlOPbqVRlOgDdbLTd3DZA29fp6Pb2+aLpaxgM4D4iepCI3lO3MgqFQvFSxuSyBBPRvQCWpPzpo8x8\nh6OYVzHzLiJaCOBeItrMzD+W01KhUCgUriBmru/Die4H8KfM/LDDe28GcJyZ/0fK3+prhEKhUHQw\nmPQTJWIAAAXdSURBVNm57F9ahuGBVGWJaAaAScx8jIhmArgGwC1p7/VpsEKhUCjCUNe22jcT0TYA\nVwC4k4jubr2+jIjubL1tCYAfE9EjAH4G4DvMfE8d+ioUCoWi5pKUQqFQKDoHTd8llQsiupaINhPR\n00T04br1iQURfYGI9hDRxsRr84joXiJ6iojuIaK5deoYAyJaSUT3tw5t/oKIPtB6vePbSETTiOhn\nRPQIET1ORB9vvd7xbUuCiCa1DtLe0fq9a9qXdlC4y9o3l4j+mYieaNnor/i2r2MJg4gmAfgMgGsB\nXADgHUR0fr1aReOLMO1J4iYA9zLzuQC+3/q9UzEM4IPM/DKYcuR/aI1Zx7eRmU8AuIqZLwFwEYCr\niOjfoQva1oY/AfA4zJZ3oLval3ZQuJva978A3MXM58PY6Gb4to+ZO/ILwJUAvpv4/SYAN9Wtl0C7\nVgPYmPh9M4DFrZ+XANhct46Cbf0WgKu7rY0AZgB4AMDLuqltAFYAuA/AVQDuaL3WTe3bAmB+22td\n0T4AcwA8l/K6V/s6NsMAsBzAtsTv21uvdRsWM/Oe1s97ACyuUxkpENFqAJfCbGjoijYSUU9rk8Ye\nAPcz8yZ0Sdta+CSADwEYS7zWTe1LOyjcLe07E8A+IvoiET1MRP+7tfvUq32dTBgvudV6NmFAx7eb\niGYB+CaAP2HmY8m/dXIbmXmMTUlqBYBfI6Kr2v7esW0jojcA2MvMG5CxFb6T29fCq5j5UgDXwZRL\nfzX5xw5v32QAlwH4HDNfBqAfbeUnl/Z1MmHsALAy8ftKmCyj27CHiJYAABEtBbC3Zn2iQERTYMji\nn5j5W62Xu6qNzHwEwJ0AXoHuadsrAdxARFsAfBXAa4jon9A97QMz72p93wfgNgCXo3vatx3AdmZ+\noPX7P8MQyG6f9nUyYTwI4BwiWk1EvQDeDuD2mnUqA7cDeFfr53fB1P07EkREAD4P4HFm/lTiTx3f\nRiJaYHeYENF0AK8DsAFd0DYAYOaPMvNKZj4TwI0AfsDMv4MuaR8RzSCi2a2f7UHhjeiS9jHzbgDb\niOjc1ktXA9gE4A54tK+jz2EQ0XUAPgVgEoDPM/PHa1YpCkT0VQCvBrAApp74FwC+DeAbAFYB2Arg\nbcx8uC4dY9DaNfQjAI9hPPX9CICfo8PbSEQXAvgyTBDWA5NBfYKI5qHD29YOIno1zJU+N3RL+4jo\nTJisAjDlm68w88e7pX0AQEQXA/hHAL0AngXwezC+07l9HU0YCoVCoagOnVySUigUCkWFUMJQKBQK\nhROUMBQKhULhBCUMhUKhUDhBCUOhUCgUTlDCUCgUCoUTlDAUCoVC4QQlDIUiASKa33oewgYi2kVE\n21s/HyOiz5TweV8ioueI6A9z3vPvWs8v2Jj1HoWiCujBPYUiA0R0M4BjzPw/S/yML8JcFf4vBe87\nA+YxxReWpYtCUQTNMBSKfBAAEFFf4ilz64joy0T0o9ZT2t5CRH/Telrb3UQ0ufW+VxDR+tZ12d+1\nl7xlfUbrf95KRBtbT+77Ydp7FIq6oIShUIThTJgHCd0A4P/CPLXsIgCDAK5v3cr7twB+k5l/CeZp\niv/NQe6fA7imdU36b5SiuUIRiMl1K6BQdCAYwN3MPEpEvwDQw8zfa/1tI8xTE8+FeeLefeaSXkwC\nsNNB9r8C+DIRfQNAbplKoagaShgKRRiGAPPQJCIaTrw+BjOvCMAmZn6lj1Bmfh8RXQ7gegAPEdEr\nmPmglNIKRQy0JKVQ+MNlPeFJAAuJ6ArAPDiKiC4oFEy0hpl/zsw3A9gH8/Q+haIR0AxDocgHJ76n\n/Qyc+lhLZuZhIvotAJ8mojkwc+2TAB7P+QwA+GsiOgeGlO5j5sdiG6BQSEG31SoUNaK1rfY7zPzN\ngvethtl+q9tqFbVBS1IKRb04AuC/Fhzc+1WYR4Xuq0wrhSIFmmEoFAqFwgmaYSgUCoXCCUoYCoVC\noXCCEoZCoVAonKCEoVAoFAonKGEoFAqFwgn/H1LdQMHLC0ZzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f16e7ee2d50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "solver = ode('cvode', rhseqn, old_api=False, one_step_compute=True)\n",
    "time = 0.\n",
    "solver.init_step(time, initx)\n",
    "plott = []\n",
    "plotx = []\n",
    "while True:\n",
    "    solution = solver.step(60)\n",
    "    if solution.errors.t:\n",
    "        print ('Error: ', solution.message, 'Error at time', solution.errors.t)\n",
    "        break\n",
    "    #we store output for plotting\n",
    "    plott.append(solution.values.t)\n",
    "    plotx.append(solution.values.y[0]) \n",
    "    if solution.values.t >= 60:\n",
    "        #back up to 60\n",
    "        solver.set_options(one_step_compute=False)\n",
    "        solution = solver.step(60)\n",
    "        plott[-1] = solution.values.t\n",
    "        plotx[-1] = solution.values.y[0]\n",
    "        break\n",
    "plt.plot(plott,plotx)\n",
    "plt.xlabel('Time [s]')\n",
    "plt.ylabel('Position [m]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By inspection of the returned times you can see how efficient the solver can solve this problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "plott length: 1328 , last computation times: [59.2297953049153, 59.28543421477497, 59.34107312463465, 59.39671203449432, 59.452350944353995, 59.50798985421367, 59.56362876407334, 59.61926767393302, 59.67490658379269, 59.730545493652365, 59.78618440351204, 59.84182331337171, 59.897462223231386, 59.95310113309106, 60.0]\n"
     ]
    }
   ],
   "source": [
    "print ('plott length:', len(plott), ', last computation times:', plott[-15:]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
